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multidimensional description is deemed necesssary. Thus, anode spotting may be 
precipitated by the elevation of T, among other factors. A review of transpiration cooling 
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I. INTRODUCTION 


Several types of space flight propulsion systems have been developed over the 
years. These include chemical, nuclear, electric and solar propulsion. The majority 
of space thrusters to date have been chemical rockets. Although the Chinese used 
rockets over 800 years ago, true development of rocket propulsion took place during 
this century [Ref. 1]. Chemical thrusters give high thrust-to-weight ratios, larger than 
unity, and have been fully developed in the form of space launch vehicles and 
attitude control thrusters. In contrast, other propulsion systems have been developed 
only to the proof-of-concept stage, and essentially remain at this stage of 
development. Nuclear propulsion was studied with the NERVA (Nuclear Engine for 
Rocket Vehicle Application) thruster in the 1960’s, and abandoned [Ref. 2:pp. 517- 
519]. Electric propulsion flights during the 1960’s included the U.S. SERT-1 (Space 
Electric Rocket Test) in 1964 and the U.S.S.R. Yantari-1 rocket in 1966. Solar- 
electric propulsion was demonstrated via the SERT-2 rocket in 1970, powering the 
electric thruster from power generated by solar cells. Further electric propulsion 
research flights in the 1980’s included the U.S. Navy’s NOVA-1 satellite in 1981, and 
Japan’s MS-T4 satellite,launched from the Space Shuttle. Beyond this, nonchemical 


thrusters have only been used in auxiliary roles, such as station-keeping and attitude 


control on geosynchronous satellites. NASA’s Project PATHFINDER in the mid- 
1980’s proposed the use of a megawatt-level electric plasma thruster for a manned 
Mars mission. However, development of this project was never funded. 

In comparing the different propulsion schemes, a primary performance indicator 
is specific impulse, defined as the ratio of thrust to the rate of propellant usage, or 
alternately, propellant effective exhaust velocity (u,), divided by the sea-level 


gravitational constant, (g,). 


ee ce Us) 





Chemical rockets are inherently limited in performance by the total energy available 
in the fuel/oxidizer combustion process, so that the total enthalpy available for 
conversion into exhaust kinetic energy is limited. Exhaust velocity is also limited by 
material heating limitations of the combustion chamber and nozzle throat, and 
"frozen flow Losses” (unrecoverable energy deposition in internal modes of the gas) 
(Ref. 3:pp. 4-5]. Peak specific impulse for liquid chemical propellants is presently on 
the order of 450 seconds. This capability is completely sufficient for the tasks of 
launch to low earth orbit (LEO), attitude control, station keeping, and such missions. 
However, for missions such as manned interplanetary exploration, chemical 


propulsion can be shown to be clearly inadequate. A comparative analysis of a Mars 


i) 


mission using chemical and electric propulsion systems shows the large difference in 
mass payload ratio (final mass/initial launch mass) for the two systems. A chemical 
system using a high impulse Hohman ellipse trajectory delivers a maximum of 
approximately 10% to 18% of the launch mass to the Martian surface [Ref. 4:p. 115]. 
In comparison, an electric system using a low impulse spiral trajectory could deliver 
from 20% to 60% of the launch mass, depending on the desired transit time. Each 
mission assumes transit from low Earth orbit to Mars orbit. An electric propulsion 
system would still need a high thrust propulsion system to reach the Martian surface 
Ref. 5:pp. 344-346]. The large difference in payload ratio is due to the much larger 
exhaust velocity and more efficient use of fuel by electric propulsion. Thus, some 
form of electric or hybrid electric thruster would seem to be in order for such 
interplanetary missions. However, due to the low thrust-to-weight ratio of electric 
thrusters, they must be launched into orbit by other means. Their usefulness is 
restricted to space thrusters, not to launch systems. 

With specific impulses of as high as 10,000 seconds, electric propulsion offers 
the performance envelope needed for manned interplanetary missions. Electric 
propulsion is divided into three types of thrusters: electrostatic, electrothermal, and 
electromagnetic. The type relevant to this work is the magnetoplasmadynamic (MPD) 
thruster, an electromagnetic propulsion system that utilizes the Lorentz force created 
by an electric current together with its induced magnetic field to propel a gas that 
has been heated to the plasma state. According to electromagnetic theory, a 


conductor carrying a current produces an induced magnetic force perpendicular to 


the current. The applied electric field and its induced magnetic field interact to 
produce the Lorentz force (Ff = 7xB) perpendicular to both fields on the 
conductors. This briefly summarizes the concept behind the "self-field" MPD 
accelerator [Ref. 2:pp. 485-486]. MPD performance is enhanced by adding magnetic 
coils to the thruster, thus strengthening the magnetic field and, as a consequence, the 
Lorentz force and thrust. This thruster is appropriately called an "applied-field" MPD 
thruster. MPD thrusters have shown specific impulses of up to 7,000 seconds and 
efficiencies as high as 70% [Ref. 6:pp. 2-3]. Performance of MPD thrusters is limited 
by several factors, including electrode erosion, current spotting, frozen flow losses, 
and electrode power deposition. Specifically, anode power deposition is the single 
largest power loss mechanism in MPD thrusters operating at submegawatt power 
levels [Ref. 7]. In the following work, we review and analytically model the MPD 


anode, including the sheath and anode potential drop. 


Il. LITERATURE REVIEW 


Anode losses significantly limit magnetoplasmadynamic (MPD) thruster 
performance. Much effort has been placed on characterizing these losses and on the 
nature of power deposition in the anode [Refs. 8-14]. As much as 80% of thruster 
total power may end up being deposited in the anode at sub-megawatt power levels 
[Refs. 8,15]. This power deposition together with current constriction at the anode 
surface present serious problems to thruster cooling and performance, as well as to 
anode lifetime. Before any practical design can be achieved, a more thorough 
understanding of the phenomena at the anode, particularly the anode sheath, must 
be gained. Studies have shown that the anode power fraction depends on thruster 
power, current, mass flow rate, and the parameter J*/m (Refs. 8,12,13,16]. It has also 
been shown that the anode fall voltage is inversely proportional to anode current 
density [Refs. 13,16]. It is believed that a better understanding of the role of an 
elevated electron temperature, of current flow dimensionality, and of current 
unsteadiness are prerequisites for the evolution of any practical MPD thruster design. 

Computer codes that accurately describe observed data from steady-state MPD 
thrusters have been developed [Refs. 17-19]. However, these codes do not adequately 
describe observed data from quasi-steady thruster experiments. It has been suggested 


that the lack of proper electrode modelling (1.e., sheaths and fall potentials) in these 


codes may explain this discrepancy [Ref. 6]. Limited analytical work has been done 
in modelling the sheath and ambipolar regions at the anode, influenced perhaps by 
the difficult set of coupled, nonlinear partial differential equations involved. Hugel 
[Ref. 12] and Subramaniam [Ref. 20] address the influence of the sheath region, but 
do not model the electric field, temperature, or sheath fall voltage. 

Given the minuscule extent of the sheath versus thruster anode curvature, the 
problem at first appears one-dimensional in nature. A one-dimensional, collisional, 
equilibrium solution can satisfactorily reproduce the observed electric field and 
charge density distributions for the entire sheath and ambipolar regions for a sheath 
where the electron temperature equals that of the heavy species [Ref. 21]. However, 
this model cannot describe any decrease in current density away from the surface, or 
current constriction, at the anode surface which might be necessary in 
nonequilibrium. A two-dimensional model, developed by Biblarz and Dolson [Ref. 
14], represents these phenomena and predicts the voltage drop in the region. It is 
shown that the sheath must account for a majority of the anode voltage drop, and 
that the sheath extent must be greater than the Debye length [Rets. 14,21]. Thus, 
a combination of one- and two-dimensional approaches appears to better describe 
sheath behavior. Incorporation of modelling of this sort may improve the ability of 
the computer codes cited above to properly describe quasi-steady thrusters. 

Next, a description of the anode region is presented in order to delineate some 


of the possible effects of temperature. 


lil, ANODE DESCRIPTION 


A. THRUSTER GEOMETRY DESCRIPTION 

The majority of plasma thrusters to date have consisted of a central cathode 
rod surrounded by an annular shell anode, as shown in Figure 1 [Ref. 23]. The 
thruster illustrated is sufficient to produce needed thrust at current levels above one 
kiloamp. Below this level, an external magnetic field produced by an annular magnet 
is needed to ensure sufficient Lorentz force on the plasma propellant to meet thrust 
requirements. [Ref. 8]. As illustrated in Figure 1, the 7xB body force simplifies into 
an axial (7,B,) body force, which provides direct electromagnetic thrust ("blowing"), 
and aradial (-3,B,) body force, which provides electromagnetic compression of the 
plasma and a subsequent pressure force along the cathode surface ("pumping"). [Ref. 
6] 

A notable exception to this geometry is the Stationary Plasma Thruster (SPT), 
a design from the former Soviet Union. The SPT is an example of a plasma 
propulsion system known as a Hall Current Plasma accelerator. An electric field 1s 
applied axially to a stream of flowing plasma, in addition to a magnetic field with a 
strong radial component, which is applied by an external electromagnet. When the 
axial electric field is applied and a current flows through the plasma, an azimuthal 


component of current is induced, Le., the "Hall" current. 
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Figure | - Magnetoplasmadynamic (MPD) Thruster, with Axial and Radial Forces 
on Plasma Indicated. [Ref. 23] 


Thrust is produced by electrostatically accelerating the ions in the plasma, as well as 
through the induced Lorentz force mentioned above. A strong radial magnetic field 
is applied to the plasma, whose properties are controlled to make the electron 
Larmor radius’ small compared to the mean free path’, while the ion Larmor radius 
is comparatively large. As a consequence, electron mobility in the axial direction is 
greatly reduced. Thus, the electric field energy is given mainly to the ions, producing 
axial ion acceleration. Collisions with neutral particles serve to accelerate the entire 
neutral plasma. [Ref. 24] 

A pair of the final prototype design developed, the SPT-100, have been 
acquired by NASA recently from Fakel Enterprise in Kaliningrad, Russia, and are 
undergoing performance evaluation at the Jet Propulsion Laboratory. Designed at 
the Kurchatov Institute of Atomic Energy (IAE) in Moscow, USSR in the 1960’s, 
smaller versions of the SPT-100 (SPT-50 and SPT-70) were flown beginning in 1972’. 
A specific impulse of 1,600 seconds and 50% efficiency, as well as space flights of 
fifty similar thrusters is claimed. The specific operational characteristics of the 
thruster are not well understood presently. Bohm diffusion of electrons and a 
phenomenon called "near-wall conductivity" have been proposed to explain the 


thruster’s operation. This thruster is shown in Figure 2. [Ref. 25] 

' Larmor radius is the radius of the helix traversed by a charged particle moving 
in a magnetic field. 

* Mean free path is the distance traveled by a particle before making a collision. 


“The suffix (i.e.,"-70") indicates the characteristic diameter of the thruster in 
millimeters. 
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Figure 2 - Stationary Plasma Thruster [Ref. 25] 


B. ELEMENTARY SHEATH FORMULAE DESCRIPTION 


1. Discussion 
Voltage losses and anode power deposition account for most of the 
inefficiency of plasma thrusters. In order to understand these losses, the anode region 
must be understood and related phenomena explained and modelled. As shown in 
Figure 3, a substantial drop in voltage occurs in a short distance from the anode 


surface. 
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Figure 3 - Electric Field Between Two Electrodes, Including "Positive Column". [Ref. 
27 | 


The anode fall region may be divided into two parts, the sheath and ambipolar 
regions. The plasma attempts to adjust itself near electrodes so as to shield the main 


body of the plasma from the electric field [Ref. 26]. The sheath is the region closest 


to the anode surface within which the ion and electron number densities are unequal, 
with the electrons dominating the region. A high electron space charge exists at the 
anode surface. This 1s caused by the anode collecting incoming electrons in 
completing the arc current with the cathode. Positive ions are produced within the 
sheath by electron impact of neutral gas molecules, and the ions are repelled toward 
the cathode. At the cathode end of the anode drop region, the density of positive 
ions is high enough to almost neutralize the electron space charge, thus forming the 
positive column or core plasma. The essential positive ion current is created in this 
way near the anode. A more complete description of this process may be found in 
Cobine [Ref. 27] and von Engel [Ref. 28]. A fundamental characteristic of plasma 
behavior is its tendency toward electrical neutrality. Whenever local charge 
concentrations arise or external potentials are introduced into a system, these are 
shielded out in a distance known as the "Debye length". This distance must be much 


smaller than the system dimension for the ionized gas to be considered a plasma 


(Ref. 29]. Equation (2) gives the Debye length [Ref. 26]. 


é kT 
a aS = 16910 s (m) (2) 
2 
ne n, 





"The Debye length effectively describes the radius of a shell around a charged 
particle outside of which the potential of the particle is not seen. 


| 


Another distance of interest is the electron mean free path, or distance traveled by 
a particle before making a collision. Equation (3) is from a derivation of Lin, Resler, 
and Kantrowitz [Refs. 30,31] giving the mean free path, with 4, being the 


approximate sheath length. 


Ce 7 a (3) 
n.(e?/3kT)'In(A.e7/3kT) 


Since the sheath extends at most a few mean-free lengths from the anode surface, 
curvature of the anode does not affect the governing equations in high pressure 
discharges. Thus, the region may be described in one dimension, the distance "y’ 
from the anode surface. While the Debye length is sometimes assumed as the sheath 
extent, Reference 22 showed that the sheath thickness is a function of the anode fall 


voltage and the electron temperature. Equation (4) gives the appropriate form. 
— €, Q, 7 do e®, (4) 
. en_ kT. 


An example case with a fall voltage of 100 volts gives a sheath extent of 





A, = 2.352x10° m. This compares to a computed Debye length of 


hp 


1.690x10°° m. Therefore, the sheath can be an order of magnitude larger 


than the Debye length. 
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Nasser [Ref. 32] discusses an elementary theoretical approach to the glow 
discharge problem. He suggests a set of four one-dimensional! ordinary differential 
equations, including the electron and ion current and number density equations, in 
addition to Poisson’s equation. Most solution attempts have failed, with the boundary 
conditions being identified as the culprit. A similar attempt for the plasma thruster 


is discussed below. 


2. Simplified Formulation 

The steady probe equations are first written [Ref. 21] in their simplest 
form. The anode is assumed to operate as a heavily biased probe, which is true for 
low enough currents when the anode is not a source of ions. Whenever the 
temperature can be considered fixed, the energy equations are implicitly satisfied 
and, since ion inertia is neglected, the resulting set consists only of two species 
continuity equations and Poisson’s equation. These equations are written in terms 
of y, which is the coordinate outward from the planar positive surface. Constants and 


variables are listed in Table 1. 
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TABLE 1 - NOMENCLATURE 


a...characteristic length of plasma n,...species number density at core plasma 
D, ,...species diffusion coefficient _N...total number density 

e...elementary charge constant T...temperature 

Fewvelectric field T,---neutral species temperature 


E,...electric field at anode surface §_@...two-body recombination coefficient 


E,...electric field at core plasma €,...permittivity constant 
Ji,e--Species current density v ...lonization coefficient 

J...total current HL; ,..Species mobility coefficient 
k...Boitzmann’s constant ®_...anode fall potential 
K...current parameter 4...mean-free distance 

N,; .--Species number density A44--Debye length 

n,...time rate-of-change of n, 4,...oheath thickness 


Note: Species subscripts denote ions (1) and electrons (e). 


I 





d 
j, = ew nE - (eD)— (5) 
dy 
d 
j, = eum E + (eD)— (6) 
dy 
eS (n - 1) (7) 
dy a 
SS Naeaill (8) 
eD 
= us 9 
ae kan ( ) 


Here, the j’s are species contributions to the total current density. The existence of 


negative charges as free electrons is pivotal in the formulation. Next, the Einstein 
relation, equation (9), is introduced to write the mobilities in terms of the diffusion 
coefficients. We assume that the diffusion coefficients remain constant in the 
problem. 

Equations (5) and (6) are next solved for dn, ./dy. The species current 
density equations are found from the net reaction rate of the plasma. Equations (10) 


and (11) combine to produce space derivatives for species current density. 
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fe = Vso, (10) 


eee ee (11) 
dy dy : 





Combining equations (5)-(11) produces a set of five coupled, non-linear differential 
equations describing the sheath. These are nondimensionalized to adjust all variables 
to the first order, and are rewritten below as equations (12)-(16), with 
nondimensionalized variables denoted by "X". Nondimensionalization can be 
accomplished as follows: The species number densities n,n, are divided by their 
values at infinity to produce output from the anode surface to unity at the ambipolar 
boundary. The current densities j,,j, are divided by the total current, allowing the 
Output to show the "mirror behavior" of the two currents. The electric field is divided 
by the initial anode value to give output starting from unity at the surface and 
decreasing to the final core field value. The variable "y" is divided by the 
characteristic length’ of the plasma, "a", producing y. These corrections allow all 


Output to vary in the range from zero to one, as a function of distance from the 


anode. 


°The characteristic length is defined so as to cancel the multiplying factor in the 
electric field equation, (14), (a= 1.107 x 10°). This allows a physical interpretation of 
the ion/electron number densities, as well as the decay rate of the electric field. 
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Attempts to solve this equation set using the computer code discussed below shows 
the set to be extremely sensitive to initial conditions. The computer code solver uses 
a "marching" scheme from the anode to the undisturbed plasma. The initial 
conditions are chosen to produce the electric field potential drop observed in actual 
thrusters. First and second space derivatives of the electric field are used as 
diagnostic checks to ensure reasonable output values and indicate instability of the 
integration process. Figure 4 shows the required resulting curves for the electric field 


and its first and second derivatives. 
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Figure 4 - Electric Field and First Two Space Derivatives Used as Diagnostic Checks 
for Integrator Output. (Plotted for a Generic Exponential Function). 
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Ecker characterizes the plasma at the anode as a double sheath, with the inner 
section called the "inertia sheath", and an outer section called the “energy loss 
section”. The inner section shows a potential rise of the order of one volt, with the 
outer section showing the exponential potential drop shown in Figure 4. While this 
double sheath may in fact describe the actual sheath region, the formulation above 
only models the potential drop portion of the sheath, and does not attempt to 
produce the potential rise of the inner sheath. In addition, Ecker’s current 
constrictions are of a "macroscopic" nature, whereas those of Reference 14 and this 
work are "microscopic" [Ref. 33]. 

Data for a 6,000°K Nitrogen plasma were used to test the equation set [Ref. 
21]. Producing a proper solution required adjusting the initial conditions to force the 


curve shapes discussed above. Using Equations (2), (3), and (4), the mean free path, 


Debye length, and approximate sheath extent are calculated as A = 9.352x107>° m, 
Ap = 1.690x10° m,and A, = 2.352x10~° m (this assumes a drop voltage of 100 


volts). 
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3. Approximate Formulation 
Reference 21 explores the above equation set by taking advantage of the 
symmetry among the equations, and introducing two parameters, K° and K~, shown 


below. 








ee (17) 
eD. eD. 

= (18) 
eD eD. 


It can be shown that the resulting equations can be manipulated to yield a single, 
ordinary differential equation for the K’s in terms of the electric field. The resulting 
equation can be scrutinized for two distinct temperature regimes. Note that while 
the total current density, J, is constant in a steady, one-dimensional case, the K’s can 


vary and will in turn also depend on the degree of reactivity of the plasma (n,), ie., 


a) (19) 
dy dy . 





Because ion diffusion is much slower than electron diffusion, it can be shown that the 


K’s are related by 


(20) 





As will be evident, at the electrode surface the K’s are equal to each other and at the 


undisturbed plasma, K~ = 0. The total current density may be evaluated from 


Zl 


J =en.v (2a 


ema 


where V,,, Is the electron drift velocity beyond the ambipolar region which is strictly 
a function of E,/N, (1.e., of the ratio of undisturbed electric field to the total number 


density). 


a. Effects of Temperature on Anode Constriction 

It is useful to investigate the overall effects of temperature. Since 
temperature will be considered constant, it comes in as a parameter in this 
formulation whereas charge density and electric field remain as variables. Intuitive 
arguments will be introduced which suggest that the electron and ion/neutral 
temperatures play a rather singular role in determining the intrinsic dimensionality 
of the problem, (i.e., there are cases when the geometry of the current lines is not 
necessarily impressed by the electrode geometry). Since the problem is described by 
moderate pressure, largely collisional sheaths, the ion and neutral temperatures are 
anticipated to remain reasonably equal. Depending on the gas, the electron 
temperature, on the other hand, can be elevated from the gas temperature at the 
anode where actual magnitudes depend on the local value of E/N. In order to get 
a perspective on the effects of temperature, we shall consider two extremes, namely, 
the case where the electron and ion temperature are the same (the equilibrium case) 
and the case where the electron temperature is substantially elevated from that of 


the ions/neutrals (the two-temperature case). 


a 


(1) Casel: T, = T, = T, (Equilibnum) 
The charge densities can be eliminated by combining equations 


(5)-(9), (17) and (18). The resulting equation can be shown to be 
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If the electric field decreases monotonically from the wall to the undisturbed plasma 
feeeronr,> E.),thenasy*>“,E-E, E-0, E”’-0. 


So that in equation (22) above the "outer solution” becomes: 

Ko Sie (23) 
Now this represents an acceptable solution from a physical point of view. Moreover, 
as iy => © ; 
n = D(K’)’ = 0 (24) 
which is also acceptable for an equilibrium situation at the undisturbed plasma. 


Results [Ref. 21] are shown in Figure 5 for the case of nitrogen at 6000°K using an 


approximate electric field distribution. 
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Figure 5 - Electric Field and Species as a Function of ¥, Distance From Anode. An 
Approximation Using a Shaped Electric Field and Isothermal Plasma [Ref. 21]. 
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(2) Case ll: T, >> T, = T, (Two-Temperature) 
In this case the same procedure as before yields the following 
equation where terms divided by T, have been dropped when compared to their 


counterparts divided by T,,. 


: 2 
Rep. 2) , “fo ppv, Sopmp: _ Sop (25) 
E kTD. KT. eE e 





(K*)’ - 


Assuming the same monotonic decrease as before for the electric field from the wall 
to the plasma proper, asy>7o,E-E,, E’-0, E”’-0. 


Then the outer solution becomes 


. 4) 
eS ee ihe (26) 
dy kT eD, 


Or, K* + (constant) y + constant, and n,, keeps increasing with y. 

This is not the proper outer solution for the one-dimensional, equilibrium 
plasma that we seek because the net ionization rate continues to increase well inside 
the plasma proper where conditions should saturate, yielding a constant electric field. 
Therefore, as formulated, Case II is not amenable to a one-dimensional solution. 
References 14 and 21 show how this case can be analyzed under a multidimensional 
approach. These references also discuss a method for describing the electron 
temperature as a function of E/N, then how to couple a simplified energy relation 
which satisfactorily describes a two-temperature plasma. The necessary ingredient 


to make equation (26) approach zero beyond the decrease of E to E, Is to allow J 


a» 


to fan out as indicated in Figure 6. Thus, in equation (26), the product "EJ" can 
bring down the charge production rate to arbitrarily low values. Alternatively, it is 


possible to explore techniques of bringing the electron temperature down to be in 


means. 


U (x,y) 
= 


closer equilibrium with the ions and neutrals. Transpiration cooling is one such 
B 


sexy} 
undisturbed plasma 


UW 


wall node 





Figure 6 - Two-Dimensional Model of Current Paths Showing Periodic Structure. 
Thermal Instabilities and Inhomogeneities Would Favor One Site Over Others and 
a Single Macroscopic Constricttlon May Then Be Produced. [Refs. 14, 34] 
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b.  Sumuilanty to Vacuum Arc Phenomena 

Instability phenomena observed in vacuum arcs [Ref. 35] are very 
similar to those observed in self-field thrusters [Ref. 12]. After the establishment of 
the current, the anode region operates in a vapor that issues from the electrodes. 
In vacuum arcs, Miller characterizes the anode region as operating in one of five 
distinct modes, ranging from a passive, low current mode to a high current, fully 
developed spot mode [Ref. 36]. Given the similarities mentioned above, vacuum arc 
anode research should be helpful in the understanding of MPD thruster transition 
to the anode spot mode. Existence diagrams after Miller [Ref. 36] are shown in 
Figure 7, which divide operating modes into regions as a function of anode current 

versus electrode geometry. Figure 7 shows the transition from glow to spot mode. 
Anode spot formation at high currents is clearly a factor in limiting anode 
lifetime. Various phenomena have been related to anode spotting. Hugel [Ref. 12] 


relates the transition to spotting mode to an increase in J*/m above a critical level. 


A separate factor connected with the spot mode 1s surface temperature of the anode. 
Rich, et.al., [Ref. 37] show that anode spotting is preceded by a luminous "footpoint’ 
and followed by local melting prior to spot formation. Separately, Schuocker [Ref. 
38] finds a connection between spotting initiation and the factors of anode 
evaporation and magnetic constriction in vacuum arcs with high currents. 
Experimental investigations must be performed to see if the above-mentioned 


vacuum arc criteria apply to self-field thrusters. 
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14 After Kaltenecker (1982) 
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Figure 7 - Anode Discharge Modes as a Function of Current and Gap Length. [Ref. 
36] 


C. COMPUTER CODE 

Rather than using linear approximations to equations (12)-(16), the nonlinear 
set was used, with initial conditions adjusted in an attempt to produce observed 
electric fields from probe data. First and second space derivatives of the electric field 
were used as diagnostic checks to ensure computed output was reasonable. Initial 
conditions computed from the approximate formulae in Reference 21 were used. The 
equation set above presents a difficult problem for two reasons, nonlinearity and 
multiple time constants. The species number density equations, (12) and (13), both 
contain a nonlinear term, each with a time constant of its own. In addition, the 
electric field equation, (14), adds a possible third time constant. This constitutes a 
"stiff’ set of equations. Attempts were made to solve the set with the data discussed 
above, using Gear’s method of backward differentiation, in hopes that the variables 
would change slowly enough with each iteration to render a convergent iterative 
process. As described in Reference 39, if some reactions are slow and others fast 
among a set of coupled equations, the fast ones will control the stability of the 
method. This is addressed in the DGEAR program available from the International 
Mathematical & Statistical Library (IMSL). The latter software contains an Adams 
predictor-corrector method, as well as Gear’s method, which is well known for its 
success at solving stiff equation sets. The DGEAR software allows for a choice of 
functional or chord iteration methods, as well as a choice of Jacobian matrices. A 
more detailed discussion of this software can be found in Reference 39 and in the 


IMSL library. [Ref. 39] 
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D. COMPUTATIONAL RESULTS 


Numerous computer runs were completed using the initial conditions taken 
from Reference 21. In addition, data for the ionization coefficient v, Figure 8, was 


taken from References 40 and 41. 
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Figure 8. fonization Coefficient v as a Function of E/N for Nitrogen for Various 
Vibrational Temperatures (Refs. 40,41). 


Various combinations of initial conditions and ionization coefficient were used. As 


mentioned above, the electric field and tts space derivatives were used as 
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diagnostic/reasonability checks on the output. Individual, as well as multiple 
computer runs were attempted to model the sheath region. Nonlinearities in the 
equation set are clearly seen in Figure 9. The ion number density does not reach that 
of the electrons, and the latter population growth rate continues to grow without 
bound. The shape of the electron population curve Is very sensitive to its initial value. 
As shown in Figure 9, the latter population has too high a growth rate when 
compared to the ion population, and the latter does not "catch up". Increasing the 
initial value of fi, flattens out this curve to a reasonable shape. Above an initial 
value of approximately 0.06, however, the plot of fi, "dips" after a certain distance 
and then continues to increase as expected. This gives an approximate upper value 
for this initial value. To avoid instabilities like this, small "slices" were taken of the 
output after a small number of integration steps and multiple runs were used to form 
a "cut and paste” plot of the region. When a reasonable plot shape was produced, the 
value of ionization coefficient was varied in the "slices" to attempt to produce the 
required end values for electric field and species population. Both multiple and 
equilibrium values for the ionization coefficient were used. When the data showed 
signs of instability and failure to follow the required forms of Figure 4, a "slice" was 
made in the data stream, and the data points from this point used to start a new 
computer run. This approach was taken in the hope of avoiding singularities in the 
integration from anode surface to ambipolar region. In addition to the diagnostic 
checks shown in Figure 4, an additional data check is provided by the transition from 


the sheath to the ambipolar region. 
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Figure 9. Species Number Density plots for Individual Computer Run, Showing 
Divergent Tracks for Ion and Electron Populations, and Effect of Nonlinearity. 


As shown in Figure 5, the species number densities are equivalent in this region, as 
are their change rate. Thus, setting Equations (12) and (15) equal to each other and 


solving for fi yields a value of 0.5 in the ambipolar region. As indicated in Figures 


ae 


10-11, the output produces the desired plot slopes for electric field and species 
number density. However, the number density plots cross long before approaching 
the required value of 0.5. In addition, neither electric field nor species number 
density approaches an equilibrium value or shows sign of levelling off. Apparently, 
the multiple time constants and nonlinear portions of the number density equations 
combine to create a seemingly intractable system. Solutions for this system may be 
possible for specific, individual initial condition sets, but the problem does not appear 
amenable to this approach in general. A one-dimensional system such as this may be 
better described through the approach of boundary layer theory or nonlinear 
dynamics and chaos. Given the effort and difficulty involved in the latter, a one- 
dimensional approach such as that modelled above does not appear useful. A 
combination of one- and two-dimensional modelling would appear to be more useful, 
as discussed in Reference 14. A one-dimensional model may be useful, but only in 
an approximation approach, with a shaped electric field, such as that used in 


Reference 21. 
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Figure !0. Electric Field as a Function of ¥, Distance From Anode, Using 
Equations (12)-(16). 
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Figure 11. Species Number Density as a Function of ¥, Distance From Anode, Using 
Equations (12)-(16). 
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E. ANODE FALL VOLTAGE 

The arc discharge operates with an anode voltage drop (V,,) which should be 
on the order of the ionization potential of the gas as it is in glow discharges. This 
voltage drop is largely non-ohmic because the anode region is typically very thin and 
space-charge controlled. Since a significant portion of the anode heating comes via 
the oncoming electrons accelerating through the fall voltage, it is of interest to 
minimize this voltage, which can range from less than ten to over 40 volts in argon. 
The question arises as to what governs this anode drop and why is it such a 
noticeable function of current? [Ref. 42] 

In collisional sheaths, the anode voltage drop is largely governed by a positive 
ion generation region which forms in front of the anode. The production of ions 
reduces the space charge density and thereby permits operation at lower voltages 
than otherwise possible. Ions are most often created by electrons within their last 
few mean-free-paths before entering the anode; these ions slowly drift away from the 
anode, thereby effectively neutralizing the space charge, at a speed proportional to 
the ratio of their mobility to that of the electrons. At the anode, the electric field 
is a maximum and the electron mean energy displays a corresponding high. 

At moderate pressures, the sheath is very thin and breakdown can be visualized 
as occurring between the undisturbed plasma (which acts as a rich source of 
electrons) and the anode surface. This arrangement has the attributes that represent 
"thermionic arc breakdown" [Ref. 43], a source of electrons which is independent of 


the breakdown itself and relatively small spacings. For gases that allow cumulative 
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ionization (with some help from the tail of the Maxwellian distribution of electron 
energies), what results is a breakdown voltage appreciably below the ionization 
potential. This then could be an explanation for the low voltage breakdown 
observations [Ref. 42]. Clearly, gases with low ionization potentials and lots of 
atomic electron energy levels are preferred (such as cesium and barium) but low- 
voltage breakdown has been observed with most gases. 

The increase of the anode fall voltage above the ionization potential has been 
related to the electron Hall parameter, since a reduction of this parameter decreases 
that voltage and corresponding losses. Control of the local magnetic field through 
the use of and array of permanent magnets as well as implementation of 
transpiration cooling (which increases the electron collision frequency) have both 
yielded some encouraging results. Because the anode fall also scales up with J*/m, 
it is conceivable that current inhomogeneities and plasma instabilities which are 
reflected in this parameter are in the picture as well. [Ref. 44] 

In summary, any possible reduction of the anode fall voltage will hinge on a 
thorough understanding of the anode region, with its associated sheath and ambipolar 
regions, where electron temperature effects, ionization effects, and magnetic field 
effects play a pivotal role. If transpiration cooling is present, then additional 
phenomena of fluid-dynamic nature may come into play. Experimental observations 
with atmospheric discharges indicate the possible presence of convective effects at 


the anode. [Ref. 45] 
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IV. TRANSPIRATION COOLING 

Transpiration cooling of the anode has often been promoted as an attractive 
means of recovering a large portion of the power deposited there. Additionally, the 
onset of melting may be minimized or even avoided by active anode cooling. Rich, 
et.al., related high anode temperatures to anode spotting [Ref. 37]. Similarly, Park 
and Choi showed that low thermal diffusion leads to erosion and, consequently, 
anode damage [Ref. 46]. Active anode cooling via transpiration is one means of 
ensuring high thermal diffusion and extending anode lifetime. Early work by 
Schoeck, et. al., [Ref. 47] showed that up to 80% of the energy deposited in the 
anode is recoverable via transpiration cooling. While this study used non-convective, 
high-intensity argon arcs, it is reasonable to assume that this effect would apply in 
MPD arcs using other propellants. Although this cooling method has not been 
studied for incorporation in plasma thrusters for some time, it has been recently 
considered as a means of cooling the fuselage of the National Aerospace Plane 
(NASP). Plasma thruster designs could undoubtedly gain from this database, and due 
consideration should be given to this cooling approach for the anode. 

For a given mass transfer flow rate, the heat flux reduction to a surface is 
inversely proportional to the molecular weight of the injected gas. Use of the 
propellant as coolant as well as fuel would eliminate the need for additional tankage 


and pumps, simplifying the design considerably. Lithium has been considered to be 
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the propellant of choice, primarily because of its low molecular weight, its favorable 
ionization potential, and its low-volume tankage properties. It has a relatively low 
first ionization potential of 5.4 eV and a high second ionization potential of 75.6 eV. 
This single ionization potential range of over 70 eV compares to approximately 20 
eV for Cesium and 27 eV for Potassium [Ref. 48]. This provides a broad temperature 
range within which only single ionization will occur. Large temperatures must be 
reached within the gas before double ionization occurs. As the gas temperature is 
increased several thousand degrees Kelvin, it undergoes ionization and disassociation. 
Thermal energy deposited can be recovered through nozzle expansion at the exit. 
However, residence time of the gas is not long enough to ensure recombination. 
Thus, the energy invested in ionization and disassociation will be lost. [Ref. 49] 
Lithium has been shown to produce specific impulse figures in excess of 7,000 
seconds at 70% efficiency in steady-state thrusters, [Ref. 50] whereas all other 
propellants have been limited to less than 3,000 seconds specific impulse at less than 
40% efficiency [Ref. 6]. Subramaniam has concluded that: 
..regenerative cooling of anodes (at the specific impulse values in the MPD 
regime) is possible only with hydrogen or with alkali-metal propellants, notably 
lithium. In the latter case, the ideal anode operating mode would be 
evaporation and ionization of the propellant on the porous or wetted anode 
surface, resulting in increased ion current fraction, reduced anode voltage fall 
and utilization of part of the anode loss energy [Ref. 51]. 
Liquid coolants, as well as reducing storage requirements, offer the advantage of 


providing latent heat of vaporization for energy disposal. However, design problems 


can occur if the liquid is allowed to vaporize within the porous structure. Problems 


39 


arise due to the abrupt increase in pressure gradient as the coolant vaporizes. Since 
coolant flow generally have three-dimensional characteristics, the flow will be 
diverted around the vapor bubbles and hot spots often develop. The technical 
practicality of using molten lithium to cool a porous tungsten anode would seem to 
be beyond current technology. On the other hand, the products of decomposition of 
hydrazine (gaseous hydrogen and ammonia) have proven to be efficient and practical 
coolants [Ref. 52]. 

Given the performance figures above, using an auxiliary coolant gas even with 
high molecular weight (e.g., NH;, N,, CH,, etc.) which could serve as a propellant 
once released from a porous hot tungsten anode surface would seem to more 
practical, vice dealing with molten lithium. Experimental studies would be needed 
to compare the approaches. Kuriki and Suzuki performed experiments with a quasi- 
steady MPD thruster to study the effect of anode gas injection (Argon). At high 
currents of up to 10 kA, increases in thrust, specific impulse, and flow discharge 
Stability were observed [Ref. 53]. 

There is some question as to the likelihood of current constriction resulting 
from anode gas injection. In such a case, swirling or circulating the propellant gas 
would help to move any footpoints that developed around the anode surface and 
prevent them from becoming fully-developed spots. Additionally, an applied 
magnetic field could serve to circulate the footpoints as well. The unique advantage 


of transpiration cooling hinges on providing effective anode cooling while supplying 
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hot propellant, but the real benefit will depend on how small the amount of coolant 
required really will be. 

Transpiration cooling has proven to be as desirable as it is challenging. It is 
complicated to implement, with associated reliability problems and difficulty of 
analytical predictions. While the production of thicker boundary layers is largely 
ineffective against the electron flux heating, the cooling itself is most efficient and a 
substantial fraction of the energy transferred to the anode is recoverable. The 
arguments of Chapter Three indicate that a reduction of the electron temperature 
in the anode would have the desirable effect of reducing the initial current spotting 
which can be conjectured to be the path that leads to anode arc spots. This electron 
temperature reduction can be done most effectively by polyatomic gases (which have 
a high 6-loss factor) emanating from the anode surface [Ref. 54]. 

The arguments relating to transpiration cooling might be summarized as 
follows: 

Favorable Outcomes 

- No separate cooling mechanism for anode required, 

- Adds "hot" propellant to exhaust "recovering" most of the electrical power loss to 
the anode, 

- Quenches T, thus likely to postpone anode spotting and reduce the heating 
associated with the electron thermal energy (SkT,/2e), 

¢ Reduces bulk convective heating, 


¢ Reduces the local electron Hall parameter by increasing the collision frequency, 
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Favorable Outcomes (cont’d.) 
- Allows for some radiation cooling from the hot tungsten surface (about 120 


watts/cm’ at 2800°K [Ref. 51]), 


- Hydrogen/ammonia gases flowing through hot porous (sintered) tungsten 


represent a compatible, proven technology. 


Unfavorable Outcomes 
- May decrease the electrical conductivity in the anode region, 


- May destabilize the ionization processes in the sheath and bring about significant 


fluctuations in the current, 

- Disrupts "cathode jet" in front of the anode with unpredictable consequences, 

- Introduces propellant which may not be hot enough, not ionized enough, or not 
in the proper place for jxB acceleration, 


- Transpiration cooling through a porous (tungsten) anode is a difficult design 


problem. 


42 


V. CONCLUSIONS AND RECOMMENDATIONS 


Plasma thrusters offer distinct advantages in terms of payload delivered for 
interplanetary missions, as well as for orbital transfer. A recent comparison 
completed by Choueiri, Kelly, and Jahn shows a mass savings of 65 tons for an 
orbital transfer from low Earth orbit to geosynchronous Earth orbit using a 
quasisteady MPD thruster as opposed to an advanced chemical thruster.’ This 
superior performance comes at the expense of low thrust-to-weight ratio and long 
transit time. However, given the large cargo/logistic requirements of a manned 
interplanetary mission, delivery of payload must be maximized. Thus, further work 
to characterize more fully thruster behavior and anode contributions in particular are 
certainly warranted. [Ref. 55] 

The "cut-and-paste" method used to generate Figures 10 and 11 is not of 
practical use as a modelling method, due to the large effort involved. It did produce 
the expected electric field and species number and current density plots near the 
anode, but failed to produce the entire sheath out to the ambipolar region. The 
nonlinearity of the equation set led to a quickly deteriorating solution. A more 
practical approach using nonlinear dynamics and/or chaos must be developed to 


model the sheath numerically. 


“This assumes a specific impulse of 2,000 seconds, 600 kW of input power, and 
a 270-day transit time. 
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Depending on the propellant mass fraction used for cooling, the transpiration 
scheme discussed above presents some rather unique advantages. A hot anode 
which uses only a small amount of propellant for cooling need not be penalized for 
any lost thrust. If in addition, we increase anode lifetime by delaying the formation 
of anode arc spots, then the scheme is all the more desirable. A decrease of the 
electron temperature in the vicinity of the anode may bring about a more 
homogeneous flow of current and a reduction in the heating effect associated with 
the high electron kinetic energy. Recovery of the heat deposited at the anode would 
be most important if the propellant fraction is high. In such case, nozzle expansion 
of the hot-propellant/coolant-gas might be implemented. 

Means of limiting anode losses through decreasing anode fall voltages were 
discussed, including the control of the local Hall parameter and the implementation 
of thermionic arc breakdown. The electrical conductivity (of a nonreacting plasma) 
could possibly decrease as a result of transpiration cooling and this might increase 


the anode fall voltage. 


Additional work needs to be done in the following areas: 
- Investigate effectiveness of nonlinear dynamics and chaos in solving sheath 
equation set, 
- Incorporate adequate one- or two-dimensional sheath modelling in quasisteady 


MPD numerical codes, 
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- Investigate the role that fluid dynamic effects play in MPD thruster anode 
discharges, 

- Investigate the effect of transpiration cooling on current and plasma stability, as 
well as on thruster performance and lifetime, 

- Determine effectiveness of transpiration cooling’s increase of the collision 
frequency parameter, 

- Compare performance of gaseous propellant/coolants versus hybrid designs with 
lithium propellant/gaseous coolant, 

- Determine if required percentage of propellant gas as coolant is practical (e.g., less 
than 10%), 

- Investigate effect of surface imperfections as focal points for current constrictions 


and as precursors to anode spotting. 
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APPENDIX A 


The following software includes the calling program, SHEATH, its two 
subroutines, FCNJ and YDOT, and the DGEAR integrator. The latter is quite 
extensive in length and includes ten subroutines, including the following: DGRST, 
DGRCS, DGRPS, DGRIN, LUDATF, LUELMF, LEQT1B, UERTST, UGETIO, and USPKD. A 
detailed discussion may be found in IMSL literature or Reference 39. 


KkKk kkk kkk KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KK KKK KKK KKK KKK KK KKK KK KR KK Kk RK KK Kk kK ke ke kK 


Program Sheath 000010 
S 000020 
C- a Calling program for DGEAR integrator. Initial conditions are 000030 
C input via READ statements and keyboard entry. Output is to data 000040 
Cc files via the DGRST subroutine. Diagnostic check of output via 000050 
S Figure 4 printed to data file from DGRST subroutine. Consult 000060 
C DGEAR comments for variable descriptions not listed below. 000061 
Ci ee 000062 
REAL E,K,EPS,T1,E FMEFINE ,DIABEANUINE, Caw = 000070 
INTEGER N, METH,MITER, INDEX, IWK(1),IER,STEP 000080 
REAL*8 X,H,Y(5) ,XEND, TOL, WK 000090 
EXTERNAL YDOT, FCNJ, DGEAR 000100 
COMMON /CONST/E,K,EPS,TI,EF,EFINF,NINF,DI,DE,VINF,C1,K1,A 000110 
e 000120 
C-----2-e-0----- Constants COO LZ: 
€ Cl and Kl are constants describing the lonization coefficient. O00t22 
G They are taken from the data plotted in Figure 8. The GO90123 
c coefficient 1s equal to the nondimensionalized electric potential 000124 
@ raised to the Kl power and then multiplied by Cl. OO@ 125 
es In this way, the lonization coefficient is allowed to vary in 000126 
Cc proportion to the strength of the electric potential. OCOlz? 
C--- eee ee eee ree O00 TZs 
WRITE (*,*)*’ Input value*fer Cl (femuaeters) =” CO0rzZe 
READ (ey see 000130 
WERE DEG) Cu O00 3a 
WRITE (*,*) “input value tor Ki termateoes)... 00012 
READ (% <1) 000133 
WRITE C*; - ro 000134 
Ec Initial conditions for species number density, electric potential 000135 
E and species current density are now input (ni,ne,E,je,ji). 000136 
Ce reece err eee OOOTse 
WRITE (*,*)*’ Input values for y(l) through) y(5) (formmaces 6s oe. 000138 
READ(*,*) y(QJ5 yit2zie ye ya 000139 
WREDE (*, “)y Cl yen y 42 y eye) 000140 
GC Following constants are for plasma described in Reference 21 000141 
G (6,000 K, Init E=20,000 V/m, Final E=1,200 V/m) 000142 
Bol Gheadg? 000143 
K= lobe Zs QO01TSC 
EPS=8 .854E-12 000160 
al=6a3 O00CT7S 
EF=2E> 000180 
EFINF=1.2E4 O001T Ie 
DI=1.724E-4 000200 
DE=1.724E-1 000210 
VIO = 2.E6 000215 
VINF = 4.93E-7 000220 
Cc 000Zs¢g 
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WN kr 


DO 


+’Error Code = 


A is plasma characteristic length which shows potential drop. 
A ((EPS*EF) /(E*NINF)) = 
A i LUJE-5 

Xx OO: 

feND = 10. 

H = le-6 

TOL = 1E-6 

Meth = 2 

MITER 
INDEX 
N=5 
mae (1) = 5 

WK = 18000. 

Pee = O 

OPEN (UNIT=8 , FILE=’ SHEATH.DAT’ , STATUS=' UNKNOWN’ ) 

CALL DGEAR2 (N, YDOT, FCNJ,X,H,Y,XEND, TOL, METH, MITER, INDEX, IWK, WK, 


1 .J07E-6 


i 
a 


+IER,STEP) 


mo 3 I=-0,N 
po =0, 100 
WREITE(*,*)J,Y (I) 
WRITE (8,1)J7,Y (TI) 
mermean(T2,FS5.1,5(5X,D9.2)) 
CONTINUE 
CONTINUE 
WRITE (*,*)’Total Steps = ’,STEP,’Final Step Size = 
LER 


f 
Pel, 


CLOSE (UNIT=8 ) 
SPEOP 
END 


Cl & Fe He ie te ee ee te ie ie te ee ee ie ke kek kk kk ke kk kk 


C DUMMY SUBROUTINE FCNJ 


CERKKEKKKKKKKEKKKKKKKKKKKKKKKKKK KKK KKK 


SUBROUTINE FCNJ (N,X,Y, PD) 
INTEGER N 

REAL Y(N) ,PD(N,N) ,X 
RETURN 

END 


CERKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKEK 


C SUBROUTINE YDOT 


CC® & oe ee ee ie i ek te ke ek ke kk kk ke kkk ke kk kK ke kk 


CYGHre? ©} 


Cer Q Q OQ 


SUBROUTINE YDOT(N,X,Y,YPRIME,eprime, eprime2) 
REAL*8 X, Y(5),YPRIME(S) ,NUI,eprime,eprime2 
REAL E,K,EPS,TI,EF,EFINF,NINF,DI,DE,VINF,C1,K1,A,B1,B2,B3,B4 
COMMON/CONST/E,K,EPS,TI,EF,EFINF,NINF,DI,DE,VIO,VINF,C1,K1,A 
w= Cl * (Y(3)**K1) 
wor = VI / VIO 


Following constants are the bracketed values in Equations 12-16. 


fers left aS a variable. 


BL = ((E*EPS)/(K*TI)) * A 

Bl= 3.86E5 * A 

B2 = ((E*EFINF) /(K*TI)) * A 

B2 = 2.32E4 * A 

B3 = ((E*NINF) /(EF*EPS)) * A 

mer =]> 9.04ES * A 

B4 = ((VINF*K*TI) /(E*DE*EFINF)) * A 

B4 = 2.62E-21 * A 

Alpha = 2-body recombination coefficient (fm. Laser Kinetics 
Handbook (AFWL-TR-74-216, 1974)) (cm3/sec) 
Alrha = 9.e-8 


4°] 


000240 


000250 
000270 
000280 
000290 
000300 
000305 
000310 
000320 
000330 
000340 
000350 
000360 
000370 
000380 
000390 
000400 
000410 
000420 
000430 
000440 
000450 
000460 
000470 
000480 
000490 
000500 
000510 


UP WN Pe 


Cn Sete = = 2 ee = 
S Ni 
YPRIME (1) =  (BU* ¥ (Gi) * 3 (30 easy 
Se Ne 
YPRIME (2) = -(B * Y(2) * ¥(3))) 434) 
Cc E 
YPRIME(3)) = Bsr (1) 2) 
C je 
YPRIME (4) = -B4 * Y¥(2) * (VIT - (ALPHA * Y(1))) 
c aa 
YPRIME(5) = Ba *°¥ (2) * (Vil ALE ee) 
G 
C---Diagnostic Check of first,second derivatives----- 
G 
eprime = y(1) - y(2) 
eprime2 = yprime(1) - yprime (2) 
Gc 
RETURN 
END 
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Cc IMSL ROUTINE NAME - DGEAR DGEA0010 
Cc DGEA0020 
C--modified to return # of steps via variable "step" in subroutine call + 

lic DGEA0040 

= ¢ COMPUTER - IBM/DOUBLE DGEA0050 

ic DGEA0060 
| C LATEST REVISION - NOVEMBER 1, 1984 DGEA0070 
ac DGEA0080 

C PURPOSE - DIFFERENTIAL EQUATION SOLVER - VARIABLE ORDER DGEAO0090 
C ADAMS PREDICTOR CORRECTOR METHOD OR DGEA0100 
Cc GEARS METHOD DGEA0110 
é DGEA0120 
C USAGE - CALL DGEAR (N, FCN, FCNJ,X,H,Y,XEND, TOL, METH, DGEA0130 
Cc MITER, INDEX, IWK, WK, IER) DGEA0140 
Cc DGEA0150 

C ARGUMENTS N - INPUT NUMBER OF FIRST-ORDER DIFFERENTIAL DGEA0160 
C EQUATIONS. DGEA0170 

nC FCN - NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS. DGEA0180 

a ¢ (INPUT) DGEA0190 
C THE SUBROUTINE ITSELF MUST ALSO BE PROVIDED DGEA0200 
€ BY THE USER AND IT SHOULD BE OF THE DGEA0210 
@ FOLLOWING FORM DGEA0220 
G SUBROUTINE FCN (N,X,Y,YPRIME) DGEA0230 
Cc REAL X,Y(N) , YPRIME (N) DGEA0240 

- Cc ; DGEA0250 

| Ze DGEA0260 
E ; DGEA0270 
CG FCN SHOULD EVALUATE YPRIME(1),...,YPRIME(N) DGEA0280 
C GIVEN N,X, AND Y(1),...,Y(N). YPRIME(I) DGEA0290 
é IS THE FIRST DERIVATIVE OF Y(I) WITH DGEA0300 
é PESPECT TO X. DGEA0 310 
e FCN MUST APPEAR IN AN EXTERNAL STATEMENT IN DGEA0320 
C THE CALLING PROGRAM AND N,X,Y(1),...,Y(N) DGEA0330 
e MUST NOT BE ALTERED BY FCN. DGEA0 340 
C FCNJ - NAME OF THE SUBROUTINE FOR COMPUTING THE DGEA0350 

- Cc JACOBIAN MATRIX OF PARTIAL DERIVATIVES. DGEA0360 

‘Be (INPUT) DGEA0370 

nm Cc THE SUBROUTINE ITSELF MUST ALSO BE PROVIDED DGEA0380 
G BY THE USER. DGEA0 390 
G IF MITER=1 IT SHOULD BE OF THE FOLLOWING DGEA0400 
C FORM DGEA0410 
G SUBROUTINE FCNJ (N,X,Y, PD) DGEA0420 
G REAL X,Y(N) , PD(N,N) DGEA0430 
G ‘ DGEA0440 
ec ; DGEA0450 
e FCNJ MUST EVALUATE PD(I,J), THE PARTIAL DGEA0460 
C DERIVATIVE OF YPRIME(I) WITH RESPECT TO DGEA0470 
C Y(J), FOR I=1,N AND J=1,N. DGEA0480 
q IF MITER= -1 IT SHOULD BE OF THE FOLLOWING DGEA0490 
Gq FORM DGEA0500 
C SUBROUTINE FCNJ (N,X,Y,PD) DGEA0510 
G REAL X,Y(N), PD(1) DGEA0520 
Gc ; DGEA05 30 
e ; DGEA0540 
C FCNJ MUST EVALUATE PD IN BAND STORAGE MODE. DGEAO0SS50 
C THAT IS, PD(N*(J-I+NLC)+1I) IS THE PARTIAL DGEA0S60 
C DERIVATIVE OF YPRIME(I) WITH RESPECT TO DGEAO0570 
c Y(J). NLC IS THE NUMBER OF LOWER DGEA0580 
© CODIAGONALS FOR THE BAND MATRIX. DGEA0590 
a FCNJ MUST APPEAR IN AN EXTERNAL STATEMENT INDGEAO600 
é THE CALLING PROGRAM AND N,X,Y(1),...,Y(N) DGEAO610 
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XEND 


TOL 


METH 


MITER 


MUST NOT BE ALTERED BY FCNJ. 

FCNJ IS USED ONLY IF MITER IS EQUAL TO 
1 OR -1. OTHERWISE A DUMMY ROUTINE CAN 
BE SUBSTITUTED. SEE REMARK 1. 


- INDEPENDENT VARIABLE. (INPUT AND OUTPUT) 


ON INPUT, X SUPPLIES THE INITIAL VALUE 
AND IS USED ONLY ON THE FIRST CALL. 
ON OUTPUT, X IS REPLACED WITH THE CURRENT 


DGEA0620 
DGEA0630 
DGEA0640 
DGEA0650 
DGEA0660 
DGEA0670 
DGEA0680 
DGEA0690 


VALUE OF THE INDEPENDENT VARIABLE AT WHICHDGEAO0700 


INTEGRATION HAS BEEN COMPLETED. 
- INPUT/OUTPUT. 
ON INPUT, H CONTAINS THE NEXT STEP SIZE iN 
A. H IS USED CNELY “OhBtHE FIRST CAnL: 
ON OUTPUT, H CONTAINS THE STEP SIZE USED 
LAST, WHETHER SUCCESSFULLY OR NOT. 


- DEPENDENT VARIABLES, VECTOR OF LENGTH N. 


(INPUT AND OUTPUT) 


ON INPUT, Y(1),..:,Y(N) SUPPLY INITIAL 
VALUES. 
ON OUTPUT, Y(1),...,Y(N) ARE REPLACED WITH 


A COMPUTED VALUE AT XEND. 


- INPUT VALUE OF X AT WHICH SOLUTION IS DESIRED 


NEXT. INTEGRATION WILL NORMALLY GO 


DGEA0710 
DGEA0720 
DGEA0730 
DGEA0740 
DGEA0750 
DGEA0760 
DGEA0770 
DGEA0780 
DGEA0790 
DGEAO800 
DGEA0810 
DGEA0820 
DGEAO 830 
DGEA0 840 


BEYOND XEND AND THE ROUTINE WILL INTERPOLATEDGEAO850 


TO X = XEND. 
NOTE THAT (X-XEND) *H MUST BE LESS THAN 
ZERO (X AND H AS SPECIFIED ON INPUT). 


- INPUT RELATIVE ERROR BOUND. TOL MUST BE 


GREATER THAN ZERO. TOL IS USED ONLY ON THE 
FIRST CALL UNLESS INDEX IS EQUAL TO -1. 
TOL SHOULD BE AT LEAST AN ORDER OF 
MAGNITUDE LARGER THAN THE UNIT ROUNDOFF 
BUT GENERALLY NOT LARGER THAN .001. 

SINGLE STEP ERROR ESTIMATES DIVIDED BY 

YMAX (I) WILL BE KEPT LESS THAN TOL IN 

ROOT-MEAN-SQUARE NORM (EUCLIDEAN NORM 

DIVIDED BY SQRT(N)). THE VECTOR YMAX OF 

WEIGHTS IS COMPUTED INTERNALLY AND STORED 

IN WORK VECTOR WK. INITIALLY YMAX(I) IS 

THE ABSOLUTE VALUE OF Y(I), WITH A DEFAULT 

VALUE OF ONE IF Y(I) IS EQUAL TO ZERO. 

THEREAFTER, YMAX(I) IS THE LARGEST VALUE 

OF THE ABSOLUTE VALUE OF Y(I) SEEN SO FAR, 

OR THE INITIAL VALUE OF YMAX(I) IF THAT IS 

LARGER. 

- INPUT BASIC METHOD INDICATOR. 

USED ONLY ON THE FIRST CALL UNLESS INDEX IS 

EQUAL TO -1. 

METH = 1, IMPLIES THAT THE ADAMS METHOD IS 
TO BE USED. 

METH = 2, IMPLIES THA@eTHE STIFF METHODS OF 
GEAR, OR THE BACKWARD DIFFERENTIATION 
FORMULAE ARE TO BE USED. 

- INPUT ITERATION METHOD INDICATOR. 

MITER = 0, IMPLIES THAT FUNCTIONAL 
ITERATION IS USED. NO PARTIAL 
DERIVATIVES ARE NEEDED. A DUMMY FCNJ 
CAN BE USED. 

MITER = 1, IMPLIES THAT THE CHORD METHOD 
IS USED WITH AN ANALYTIC JACOBIAN. FOR 
THIS METHOD, THE USER SUPPLIES 
SUBROUTINE FCNJ. 
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DGEA08 60 
DGEA0870 
DGEA08 80 
DGEA0890 
DGEA0900 
DGEA0910 
DGEA0920 
DGEA0930 
DGEA0940 
DGEAO0950 
DGEA0960 
DGEA0970 
DGEAO09 80 
DGEA0990 
DGEA1000 
DGEA1010 
DGEA1020 
DGEA1030 
DGEA1040 
DGEA1050 
DGEA1060 
DGEA1070 
DGEA1080 
DGEA1090 
DGEA1100 
DGEA1110 
DGEA1120 
DGEA1130 
DGEA1140 
DGEA1150 
DGEA1160 
DGEA1170 
DGEA1180 
DGEA1190 
DGEA1200 
DGEA1210 
DGEA1220 
DGEA1230 
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MITER = 2, IMPLIES THAT THE CHORD METHOD DGEA1240 
IS USED WITH THE JACOBIAN CALCULATED DGEA1250 
INTERNALLY BY FINITE DIFFERENCES. DGEA1260 
A DUMMY FCNJ CAN BE USED. DGEA1270 

MITER = 3, IMPLIES THAT THE CHORD METHOD DGEA12 80 
IS USED WITH THE JACOBIAN REPLACED BY DGEA1290 
A DIAGONAL APPROXIMATION BASED ON A DGEA1 300 
DIRECTIONAL DERIVATIVE. DGEA1310 
A DUMMY FCNJ CAN BE USED. DGEA1320 

PETE Ee— “ere 2, ihe LIES USE THE SAME DGEA1330 
METHOD AS FOR MITER= 1 OR 2, RESPECTIVELY , DGEA1340 
BUT USING A BANDED JACOBIAN MATRIX. IN DGEA1350 
THESE TWO CASES BANDWIDTH INFORMATION DGEA1360 
MUST BE PASSED TO DGEAR THROUGH THE DGEA1370 
COMMON BLOCK DGEA13 80 

COMMON /DBAND/ NLC,NUC DGEA1390 

WHERE NLC=NUMBER OF LOWER CODIAGONALS DGEA1400 

NUC=NUMBER OF UPPER CODIAGONALS DGEA1410 

INDEX - INPUT AND OUTPUT PARAMETER USED TO INDICATE DGEA1420 

THE TYPE OF CALL TO THE SUBROUTINE. ON DGEA1430 

OUTPUT INDEX IS RESET TO 0 IF INTEGRATION DGEA1440 

WAS SUCCESSFUL. OTHERWISE, THE VALUE OF DGEA1450 

INDEX IS UNCHANGED. DGEA1460 

ON INPUT, INDEX = 1, IMPLIES THAT THIS IS THE DGEA1470 
FIRST CALL FOR THIS PROBLEM. DGEA1480 
ON INPUT, INDEX = 0, IMPLIES THAT THIS IS NOT DGEA1490 
THE FIRST CALL FOR THIS PROBLEM. DGEA1500 
ON INPUT, INDEX = -1, IMPLIES THAT THIS IS NOTDGEA1510 
THE FIRST CALL FOR THIS PROBLEM, AND THE DGEA1520 
USER HAS RESET TOL. DGEA1530 
ON INPUT, INDEX = 2, IMPLIES THAT THIS IS NOT DGEA1540 


THE FIRST CALL FOR THIS PROBLEM. INTEGRATIONDGEAI1550 
IS TO CONTINUE AND XEND IS TO BE HIT EXACTLYDGEA1560 


(NO INTERPOLATION IS DONE). THIS VALUE OF 
INDEX ASSUMES THAT XEND IS BEYOND THE 
CURRENT VALUE OF X. 

ON INPUT, INDEX = 3, IMPLIES THAT THIS IS NOT 


DGEA1570 
DGEA1580 
DGEA1590 
DGEA1600 


THE FIRST CALL FOR THIS PROBLEM. INTEGRATIONDGEA1610 
IS TO CONTINUE AND CONTROL IS TO BE RETURNEDDGEA1620 


TO THE CALLING PROGRAM AFTER ONE STEP. XEND 
IS IGNORED. 
IWK - INTEGER WORK VECTOR OF LENGTH N. USED ONLY IF 
METER y=) OR 2 
WK - REAL WORK VECTOR OF LENGTH 4*N+NMETH+NMITER. 
THE VALUE OF NMETH DEPENDS ON THE VALUE OF 
METH . 


IF METH IS EQUAL TO 1, 
NMETH IS EQUAL TO N*13. 
IF METH IS EQUAL TO 2, 
NMETH IS EQUAL TO N*6. 
THE VALUE OF NMITER DEPENDS ON THE VALUE OF 
MITER. 
THeMOITER ESeEGUAT TO 1 OR 2, 
NMITER IS EQUAL TO N*(N+1) 
TEeMOTER IS@HOUAtT Te -1 OR -2, 
NMITER IS EQUAL TO (2*NLC+NUC+3) *N 
WHERE NLC=NUMBER OF LOWER CODIAGONALS 
NUC=NUMBER OF UPPER CODIAGONALS 
IF MITER IS EQUAL TO 3, 
NMITER IS EQUAL TO N. 
IF MITER IS EQUAL TO 0, 
NMITER IS EQUAL TO 1. 
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DGEA1630 
DGEA1640 
DGEA1650 
DGEA1660 
DGEA1670 
DGEA16 80 
DGEA1690 
DGEA1700 
DGEA1710 
DGEA1720 
DGEA1730 
DGEA1740 
DGEA1750 
DGEA1760 
DGEA1770 
DGEA1780 
DGEA1790 
DGEA1 800 
DGEA1810 
DGEA1820 
DGEA1 830 
DGEA1840 
DGEA1850 


Melo. OOOO ODAND AON NN DON OODADRAOOO000O900 00090 000000 600 OOOO 2 OO Gee 


WK MUST REMAIN UNCHANGED BETWEEN SUCCESSIVE 

CALLS DURING INTEGRATION. 

IER - ERROR PARAMETER. (OUTPUT) 
WARNING ERROR 

IER = 33, IMPLIES THAT X+H WILL EQUAL X ON 
THE NEXT STEP. THIS CONDITION DOES NOT 
FORCE THE ROUTINE TO HALT. HOWEVER, IT 
DOES INDICATE ONE OF TWO CONDITIONS. 
THE USER MIGHT BE REQUIRING TOO MUCH 
ACCURACY VIA THE INPUT PARAMETER TOL. 
IN THIS CASE THE USER SHOULD CONSIDER 
INCREASING THE VALUE OF TOL. THE OTHER 
CONDITION WHICH MIGHT GIVE RISE TO THIS 
ERROR MESSAGE IS THAT THE SYSTEM OF 
DIFFERENTIAL EQUATIONS BEING SOLVED 
IS STIFF (EITHER IN GENERAL OR OVER 
THE SUBINTERVAL OF THE PROBLEM BEING 
SOLVED AT THE TIME OF THE ERROR). IN 
THIS CASE THE USER SHOULD CONSIDER 
USING A NONZERO VALUE FOR THE INPUT 
PARAMETER MITER. 

WARNING WITH FIX ERROR 

IER = 66, IMPLIES THAT THE ERROR TEST 
FAILED. H WAS REDUCED BY .1 ONE OR MORE 
TIMES AND THE STEP WAS TRIED AGAIN 
SUCCESSFULLY. 

IER = 67, IMPLIES THAT CORRECTOR 
CONVERGENCE COULD NOT BE ACHIEVED. 

H WAS REDUCED BY .1 ONE OR MORE TIMES AND 
THE STEP WAS TRIED AGAIN SUCCESSFULLY. 
TERMINAL ERROR 

IER = 132, IMPLIES THE INTEGRATION WAS 
HALTED AFTER FAILING TO PASS THE ERROR 
TEST EVEN AFTER REDUCING H BY A FACTOR 
OF 1.0E10 FROM ITS INITIAL VALUE. 

SEE REMARKS. 

IER = 133, IMPLIES THE INTEGRATION WAS 
HALTED AFTER FAILING TO ACHIEVE 
CORRECTOR CONVERGENCE EVEN AFTER 
REDUCING H BY A FACTOR OF 1.0E10 FROM 
ITS INITIAL VALUE. SEE REMARKS. 

IER = 134, IMPLIES THAT AFTER SOME INITIAL 


DGEA1860 
DGEA1870 
DGEA1 880 
DGEA1890 
JJGEA19 00 
DGEA1910 
DGEA1920 
DGEA1930 
DGEA1940 
DGEA1950 
DGEA19 60 
DGEA1970 
DGEA19 80 
DGEA1990 
DGEA2 000 
DGEA2010 
DGEA2020 
DGEA2 030 
DGEA2 040 
DGEA2050 
DGEA2 060 
DGEA2 070 
DGEA2080 
DGEA2 090 
DGEA2100 
DGEA2110 
DGEA2120 
DGEA2130 
DGEA2140 
DGEA2150 
DGEA2160 
DGEA2170 
DGEA2180 
DGEA2190 
DGEA2 200 
DGEA2210 
DGEA2220 
DGEA2230 
DGEA2 240 
DGEA2 250 
DGEA2260 
DGEA2270 


SUCCESS, THE INTEGRATION WAS HALTED EITHERDGEA2280 


BY REPEATED ERROR TEST FAILURES OR BY 
A TEST ON TOL. SEE REMARKS. 

IER = 135, IMPLIES THAT ONE OF THE INPUT 
PARAMETERS N,X,H,XEND,TOL,METH, MITER, OR 
INDEX WAS SPECIFIED INCORRECTLY. 

IER = 136, IMPLIES THAT INDEX HAD A VALUE 
OF -1 ON INPUT, BUT THE DESIRED CHANGES 
OF PARAMETERS WERE NOT IMPLEMENTED 
BECAUSE XEND WAS NOT BEYOND X. 
INTERPOLATION TO X = XEND WAS PERFORMED. 
TO TRY AGAIN, SIMPLY CALL AGAIN WITH 
INDEX EQUAL TO -1 AND A NEW VALUE FOR 


XEND. 
Step - # of integration steps taken 
PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 


- SINGLE/H36,H48,H60 


REQD. IMSL ROUTINES - DGRCS,DGRIN,DGRPS,DGRST, LUDATF, LUELMF, LEQT1B, 


oe 


DGEA2 290 
DGEA2 300 
DGEA2 310 
DGEA2320 
DGEA2 330 
DGEA2340 
DGEA2350 
DGEA2360 
DGEA2370 
DGEA23 80 
DGEA2 390 
DGEA2 400 
DGEA2410 
+ 
DGEA2420 
DGEA2 430 
DGEA2440 
DGEA2450 
DGEA2460 
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NOTATION 


REMARKS 


Le 


VERTST, UGETIO 


- INFORMATION ON SPECIAL NOTATION AND 
CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 


THE EXTERNAL SUBROUTINE FCNJ IS USED ONLY WHEN 
INPUT PARAMETER MITER IS EQUAL TO 1 OR -1. OTHERWISE, 
A DUMMY FUNCTION CAN BE USED. THE DUMMY SUBROUTINE 
SHOULD BE OF THE FOLLOWING FORM 

SUBROUTINE FCNJ (N,X,Y,PD) 

INTEGER N 

REAL Y(N),PD(N,N) ,X 

RETURN 

END 
AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED 
(IER=0) AND A NORMAL CONTINUATION IS DESIRED, SIMPLY 
RESET XEND AND CALL DGEAR AGAIN. ALL OTHER 
PARAMETERS WILL BE READY FOR THE NEXT CALL. A CHANGE 
OF PARAMETERS WITH INDEX EQUAL TO -1 CAN BE MADE 
AFTER EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. 
THE COMMON BLOCKS /DBAND/ AND /GEAR/ NEED TO BE 
PRESERVED BETWEEN CALLS TO DGEAR. IF IT IS NECESSARY 
FOR THE COMMON BLOCKS TO EXIST IN THE CALLING PROGRAM 
THE FOLLOWING STATEMENTS SHOULD BE INCLUDED 

COMMON /DBAND/ NLC, NUC 

COMMON /GEAR/ DUMMY (48) , SDUMMY (4) , IDUMMY (38) 


WHERE DUMMY, SDUMMY, AND IDUMMY ARE VARIABLE NAMES NOT 


USED ELSEWHERE IN THE CALLING PROGRAM. (FOR DOUBLE 


DGEA2470 
DGEA24 80 
DGEA2490 
DGEA2500 
DGEA2510 
DGEA2520 
DGEA2530 
DGEA2540 
DGEA25S0 
DGEA2560 
DGEA2570 
DGEA25 80 
DGEA2590 
DGEA2 600 
DGEA2 610 
DGEA2620 
DGEA2 630 
DGEA2 640 
DGEA2650 
DGEA2660 
DGEA2670 
DGEA2680 
DGEA2690 
DGEA2700 
DGEA2710 
DGEA2720 
DGEA2730 
DGEA2740 
DGEA2750 


PRECISION DUMMY IS TYPE DOUBLE AND SDUMMY IS TYPE REAL) DGEA2760 


THE CHOICE OF VALUES FOR METH AND MITER MAY REQUIRE 
SOME EXPERIMENTATION, AND ALSO SOME CONSIDERATION OF 
THE NATURE OF THE PROBLEM AND OF STORAGE REQUIREMENTS. 
THE PRIME CONSIDERATION IS STIFFNESS. IF 
THE PROBLEM IS NOT STIFF, THE BEST CHOICE IS PROBABLY 
METH = 1 WITH MITER = 0. IF THE PROBLEM IS STIFF TOA 
SIGNIFICANT DEGREE, THEN METH SHOULD BE 2 AND MITER 
SHOULD BE 1,2,-1,-2 OR 3. 
OF THE INHERENT TIME CONSTANTS OF THE PROBLEM, WITH 
WHICH TO PREDICT ITS STIFFNESS, ONE WAY TO DETERMINE 
THIS IS TO TRY METH = 1 AND MITER = O FIRST, AND LOOK 
AT THE BEHAVIOR OF THE SOLUTION COMPUTED AND THE STEP 
SIZES USED. IF THE TYPICAL VALUES OF H ARE MUCH 
SMALLER THAN THE SOLUTION BEHAVIOR WOULD SEEM TO 
REQUIRE (THAT IS, MORE THAN 100 STEPS ARE TAKEN OVER 
AN INTERVAL IN WHICH THE SOLUTIONS CHANGE BY LESS 
THAN ONE PERCENT), THEN THE PROBLEM IS PROBABLY STIFF 
AND THE DEGREE OF STIFFNESS CAN BE ESTIMATED FROM THE 
VALUES OF H USED AND THE SMOOTHNESS OF THE SOLUTION. 
IF THE DEGREE OF STIFFNESS IS ONLY SLIGHT, IT MAY BE 
THAT METH=1 IS MORE EFFICIENT THAN METH=2. 
EXPERIMENTATION WOULD BE REQUIRED TO DETERMINE THIS. 
REGARDLESS OF METH, THE LEAST EFFECTIVE VALUE OF 
MITER IS 0, AND THE MOST EFFECTIVE IS 1,-1,2,OR -2. 
MITER = 3 IS GENERALLY SOMEWHERE IN BETWEEN. SINCE 
THE STORAGE REQUIREMENTS GO UP IN THE SAME ORDER AS 
EFFECTIVENESS, TRADE-OFF CONSIDERATIONS ARE 
NECESSARY. FOR REASONS OF ACCURACY AND SPEED, THE 
CHOICE OF ABS (MITER) =1 IS GENERALLY PREFERRED TO 
ABS (MITER) =2, UNLESS THE SYSTEM IS FAIRLY COMPLICATED 
(AND FCNJ IS THUS NOT FEASIBLE TO CODE). THE 
ACCURACY OF THE FCNJ CALCULATION CAN BE CHECKED BY 
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IF THE USER HAS NO KNOWLEDGE 


DGEA2770 
DGEA2780 
DGEA2790 
DGEA2 800 
DGEA2 810 
DGEA2 820 
DGEA2 830 
DGEA2 840 
DGEA2 850 
DGEA2 860 
DGEA2 870 
DGEA2 880 
DGEA2890 
DGEA2900 
DGEA2910 
DGEA2920 
DGEA2930 
DGEA2940 
DGEA2950 
DGEA2960 
DGEA2970 
DGEA29 80 
DGEA2990 
DGEA3 000 
DGEA3010 
DGEA3 020 
DGEA3030 
DGEA3 040 
DGEA3050 
DGEA3060 
DGEA3070 
DGEA3080 
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COMPARISON OF THE JACOBIAN WITH THAT GENERATED WITH 
ABS (MITER) = IF THE JACOBIAN MATRIX IS SIGNIFICANTLY 
DIAGONALLY DOMINANT, THEN THE OPTION MITER = 3 IS 


LIKELY TO BE NEARLY AS EFFECTIVE AS ABS (MITER) =1 OR 2, 


AND WILL SAVE CONSIDERABLE STORAGE AND RUN TIME. 

IT IS POSSIBLE, AND POTENTIALLY QUITE DESIRABLE, TO 
USE DIFFERENT VALUES OF METH AND MITER IN DIFFERENT 
SUBINTERVALS OF THE PROBLEM. FOR EXAMPLE, IF THE 
PROBLEM IS NON-STIFF INITIALLY AND STIFF LATER, 
METH = 1 AND MITER = 0 MIGHT BE SET INITIALLY, AND 
METH = 2 AND MITER = 1 LATER. 

THE INITIAL VALUE OF THE STEP SIZE, H, SHOULD BE 
CHOSEN CONSIDERABLY SMALLER THAN THE AVERAGE VALUE 
EXPECTED FOR THE PROBLEM, AS THE FIRST-ORDER METHOD 
WITH WHICH DGEAR BEGINS IS NOT GENERALLY THE MOST 
EFFICIENT ONE. HOWEVER, FOR THE FIRST STEP, AS FOR 
EVERY STEP, DGEAR TESTS FOR THE POSSIBILITY THAT 
THE STEP SIZE WAS TOO LARGE TO PASS THE ERROR TEST 
(BASED ON TOL), AND IF SO ADJUSTS THE STEP SIZE 
DOWN AUTOMATICALLY. THIS DOWNWARD ADJUSTMENT, IF 
ANY, IS NOTED BY IER HAVING THE VALUES 66 OR 67, 
AND SUBSEQUENT RUNS ON THE SAME OR SIMILAR PROBLEM 
SHOULD BE STARTED WITH AN APPROPRIATELY SMALLER 
VALUE OF H. 

SOME OF THE VALUES OF INTEREST LOCATED IN THE 
COMMON BLOCK /GEAR/ ARE 

A. HOUSED, THE STEP,SIZE H LAST USED SUCCESSFULLY 


(DUMMY (8) ) 

B. NQUSED, THE ORDER LAST USED SUCCESSFULLY 
(IDUMMY (6) ) 

C. NSTEP, THE CUMULATIVE NUMBER OF STEPS TAKEN 
( IDUMMY (7) ) 

D. NFE, THE CUMULATIVE NUMBER OF FCN EVALUATIONS 
(IDUMMY (8) ) 

E. NJE, THE CUMULATIVE NUMBER OF JACOBIAN 


EVALUATIONS, AND HENCE ALSO OF MATRIX LU 
DECOMPOSITIONS (IDUMMY (9) ) 


THE NORMAL USAGE OF DGEAR MAY BE SUMMARIZED AS FOLLOWS 


A. SET THE INITIAL VALUES IN Y. 
B. SET N, 4, H, TOL, METH, AND MITE. 


SET XEND TO THE FIRST OUTPUT POINT, AND INDEX TO 1. 


. EXIT IF IER IS GREATER THAN 128. 
. OTHERWISE, DO DESIRED OUTPUT OF Y. 
G. EXIT IF THE PROBLEM IS FINISHED. 


Ge 
D. CALL DGEAR 
E 
EF 


H. OTHERWISE, RESET XEND TO THE NEXT OUTPUT POINT, AND 


RETURN TO STEP D. 
THE ERROR WHICH IS CONTROLLED BY WAY OF THE PARAMETER 


TOL IS AN ESTIMATE OF THE LOCAL TRUNCATION ERROR, THAT 


IS, THE ERROR COMMITTED ON TAKING A SINGLE STEP WITH 


THE METHOD, STARTING WITH DATA REGARDED AS EXACT. THIS 


IS TO BE DISTINGUISHED FROM THE GLOBAL TRUNCATION 
ERROR, WHICH IS THE ERROR IN ANY GIVEN COMPUTED VALUE 
OF Y(X) AS A RESULT OF THE LOCAL TRUNCATION ERRORS 
FROM ALL STEPS TAKEN TO OBTAIN Y(X). THE LATTER ERROR 
ACCUMULATES IN A NON-TRIVIAL WAY FROM THE LOCAL 
ERRORS, AND IS NEITHER ESTIMATED NOR CONTROLLED BY 


THE ROUTINE. SINCE IT IS USUALLY THE GLOBAL ERROR THAT 


A USER WANTS TO HAVE UNDER CONTROL, SOME 
EXPERIMENTATION MAY BE NECESSARY TO GET THE RIGHT 
VALUE OF TOL TO ACHIEVE THE USERS NEEDS. IF THE 
PROBLEM IS MATHEMATICALLY STABLE, AND THE METHOD USED 
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DGEA3090 
DGEA3100 
DGEA3110 
DGEA3120 
DGEA3130 
DGEA3140 
DGEA3150 
DGEA3160 
DGEA3170 
DGEA3180 
DGEA3190 
DGEA3200 
DGEA3210 
DGEA3 220 
DGEA3230 
DGEA3 240 
DGEA3250 
DGEA3 260 
DGEA3270 
DGEA32 80 
DGEA3290 
DGEA3 300 
DGEA3310 
DGEA3320 
DGEA3 330 
DGEA3 340 
DGEA3 350 
DGEA3 360 
DGEA3370 
DGEA3 380 
DGEA3390 
DGEA3400 
DGEA3410 
DGEA3420 
DGEA3430 
DGEA3440 
DGEA3450 
DGEA3460 
DGEA3470 
DGEA34 80 
DGEA3490 
DGEA3500 
DGEA3510 
DGEA3520 
DGEA3530 
DGEA3540 
DGEA3550 
DGEA3560 
DGEA3570 
DGEA3580 
DGEA3590 
DGEA3600 
DGEA3610 
DGEA3620 
DGEA3630 
DGEA3 640 
DGEA3650 
DGEA3660 
DGEA3670 
DGEA3 680 
DGEA3690 
DGEA3700 
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SUBROUTINE DGEAR 


INTEGER 


DOUBLE PRECISION 


INTEGER 


REAL 


IS APPROPRIATELY STABLE, THEN THE GLOBAL ERROR AT A 
GIVEN X SHOULD VARY SMOOTHLY WITH TOL IN A MONOTONE 
INCREASING MANNER. 

fp toe ROUWIINE RETURNS WITH IER VALUES OF 132, 233, 
OR 134, THE USER SHOULD CHECK TO SEE IF TOO MUCH 
ACCURACY IS BEING REQUIRED. THE USER MAY WISH TO 
SET TOL TO A LARGER VALUE AND CONTINUE. ANOTHER 
POSSIBLE CAUSE OF THESE ERROR CONDITIONS IS AN 
ERROR IN THE CODING OF THE EXTERNAL FUNCTIONS FCN 
OR FCNJ. IF NO ERRORS ARE FOUND, IT MAY BE NECESSARY 
TO MONITOR INTERMEDIATE QUANTITIES GENERATED BY THE 


DGEA3710 
DGEA3720 
DGEA3730 
~GEA3 740 
DGEA3 750 
DGEA3760 
DGEA3770 
DGEA3780 
DGEA3790 
DGEA3 800 
DGEA3 810 


ROUTINE. THESE QUANTITIES ARE STORED IN THE WORK VECTORDGEA3 820 
WK AND INDEXED BY SPECIFIC ELEMENTS IN THE COMMON BLOCKDGEA3 830 


/GEAR/. IF IER IS 132 OR 134, THE COMPONENTS CAUSING 
THE ERROR TEST FAILURE CAN BE IDENTIFIED FROM LARGE 
VALUES OF THE QUANTITY 

WK (IDUMMY (11) +I) /WK(I), FOR I=1, ,N. 
ONE CAUSE OF THIS MAY BE A VERY SMALL ‘BUT NONZERO 
INITIAL VALUE OF ABS(Y(I)). 
IF IER IS 133, SEVERAL POSSIBILITIES EXIST. 
IT MAY BE INSTRUCTIVE TO TRY DIFFERENT VALUES OF MITER 
ALTERNATIVELY, THE USER MIGHT MONITOR SUCCESSIVE 
CORRECTOR ITERATES CONTAINED IN WK(IDUMMY(12)+I), FOR 
I=1,...,N. ANOTHER POSSIBILITY MIGHT BE TO MONITOR 
THE JACOBIAN MATRIX, IF ONE IS USED, STORED, BY 
COLUMN, IN WK(IDUMMY(10)+1I), FOR I=1, ,N*¥N IF 
ABS (MITER) IS EQUAL TO 1 OR 2, OR FOR I=1, oN LE 
MITER IS EQUAL TO 3. 


INC. ALL RIGHTS RESERVED. 


= 19384 BY IMS, 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


(N, FCN, FCN, XH, Y, SBND; TOL, METH MITER, INDEX, 
IWK, WK, IER, step) 
SPECIFICATIONS FOR ARGUMENTS 
N,METH,MITER, INDEX, IWK(1),IER,step 
i, a, Y (N), XEND, TOL, WK (1) 
SPECIFICATIONS FOR LOCAL VARIABLES 
NERROR, NSAVE1,NSAVE2 ,NPW,NY,NC,MFC,KFLAG, 
JSTART,NSQ,NQUSED, NSTEP,NFE,NJE,1I,N0O,NHCUT,KGO 
JER, KER, NN, NEQUIL, IDUMMY (21) ,NLC,NUC 
SDUMMY (4) 


DOUBLE PRECISION T,HH,HMIN, HMAX, EPSC,UROUND,EPSJ,HUSED,TOUTP, 


EXTERNAL 


COMMON /DBAND / 
COMMON /GEAR/ 


DATA 


IF (MITER.GE.0O) 


KER = 0 
mek = 0 


URCUND 


AYI,D,DN,SEPS, DUMMY (39) 
PEN, PCN 
NLC , NUC 
T,HH, MIN, HMAX, EPSC, UROUND, EPSJ,HUSED, DUMMY, 
TOUTP , SDUMMY,NC,MFC,KFLAG, JSTART,NSQ,NQUSED, 
NSTEP, NFE, NJE, NPW, NERROR, NSAVE1, NSAVE2,NEQUIL, 
ive , LDUMMY , NO! NHCUT 
SEPS/Z3410000000000000/ 
FIRST EXECUTABLE STATEMENT 
NLC 


oe 


SEPS 
COMPUTE WORK VECTOR INDICIES 


Sy) 


DGEA3 840 
DGEA3 850 
DGEA3 860 
DGEA3 870 
DGEA3 8 80 
DGEA3 890 
DGEA3900 


. DGEA3910 


DGEA3920 
DGEA3930 
DGEA3940 
DGEA3950 
DGEA3960 
DGEA3970 
DGEA39 80 
DGEA3990 
DGEA4000 
DGEA4010 
DGEA4020 
DGEA4030 
DGEA4 040 
DGEA4050 


DGEA4070 
DGEA4 080 
+ 
DGEA4100 
“S 
DGEA4120 
DGEA4130 
DGEA4140 


, DGEA4150 


DGEA4160 
DGEA4170 
DGEA4180 
DGEA4190 
DGEA4200 
DGEA4210 
DGEA4220 
DGEA4230 
DGEA4240 
DGEA4 250 
DGEA4260 
DGEA4270 
DGEA42 80 
DGEA4290 
DGEA4 300 
DGEA4 310 
DGEA4 320 


WOW AON ae 


10 


15 


20 


NERROR 
NSAVE1 = NERROR+N 

NSAVE2 = NSAVE1+4N 

NY = NSAVE2+N 

IF (METH.EQ.1) NEQUIL = NY+13*N 

IF (METH.EQ.2) NEQUIL = NY+6*N 

NPW = NEQUIL + N 

IF (MITER.EQ.0.OR.MITER.EQ.3) NPW = NEQUIL 
MFC = 10*METH+IABS (MITER) 


N 


CHECK FOR INCORRECT INPUT PARAMETERS 


IF (MITER.LT.-2.OR.MITER.GT.3) GO TO 85 

IF (METH.NE.1.AND.METH.NE.2) GO TO 85 

IF (TOL.LE,0.DO) GO TO 65 

IF (N.LE.0) GO TO 85 

IF ((X-XEND) *H.GE.0.D0) GO TO 85 

IF (INDEX.EQ.0) GO TO 10 

IF (INDEX.EQ.2) GO TO 15 

IF (INDEX.EQ.-1) GO TO 20 

IF (INDEX.EQ.3) GO TO 25 

IF (INDEX.NE.1) GO TO 85 

IF INITIAL VALUES OF YMAX OTHER THAN 

THOSE SET BELOW ARE DESIRED, THEY 
SHOULD BE SET HERE. ALL YMAX (TI) 
MUST BE POSITIVE. IF VALUES FOR 
HMIN OR HMAX, THE BOUNDS ON 
DABS (HH), OTHER THAN THOSE BELOW 
ARE DESIRED, THEY SHOULD BE SET 


BELOW. 
DO 5 I=1,N 
WK(I) = DABS (Y(I)) 
IF (WK(I) .EQ.0.D0) WK(I) = 1.D0 
WK (NY+I) = Y(I) 
CONTINUE 
NC = N 
= xX 
HH = H 
IF ((T+HH) .EQ.T) KER = 33 
HMIN = DABS (H) 
HMAX = DABS (X-XEND) *10.D0 
EPSC = TOL 
JSTART = 0 
NO = N 


NSQ = NO*NO 
EPSJ = DSORT (UROUND) 


NHCUOT = 0 

DUMMY (2) = 1.0D0 
DUMMY (14) = 1.0D0 
GO TO 30 


TOUTP IS THE PREVIOUS VALUE OF XEND 
FOR USE IN HMAX. 
HMAX = DABS (XEND-TOUTP) *10.DO 
GO TO 4s 


HMAX = DABS (XEND-TOUOTP) *10.D0 
IF ((T-XEND)*HH.GE.0.D0) GONiTG@ees 
GO TO 50 


IF ((T-XEND) *HH.GE.0.D0) GO TO 90 
JSTART = -1 

NC = N 

EPSC = TOL 
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DGEA4 330 
DGEA4 340 
DGEA4 350 
DGEA4 360 
DGEA4370 
DGEA43 80 
DGEA4390 
DGEA4400 
DGEA4410 
DGEA4420 
DGEA4430 
DGEA4440 
DGEA4450 
DGEA4460 
DGEA4470 
DGEA44 80 
DGEA4490 
DGEA4500 
DGEA4510 
DGEA4520 
DGEA4530 
DGEA4540 
DGEA4550 
DGEA4560 
DGEA4570 
DGEA4580 
DGEA4590 
DGEA4600 
DGEA4610 
DGEA4 620 
DGEA4 630 
DGEA4 640 
DGEA4650 
DGEA4 660 
DGEA4 670 
DGEA4680 
DGEA4 690 
DGEA4 700 
DGEA4710 
DGEA4720 
DGEA4730 
DGEA4740 
DGEA4750 
DGEA4 760 
DGEA4770 
DGEA4 780 
DGEA4 790 
DGEA4 800 
DGEA4 810 
DGEA4 820 
DGEA4 830 
DGEA4 840 
DGEA4 850 
DGEA4 860 
DGEA4 870 
DGEA4 880 
DGEA4 890 
DGEA4900 
DGEA4910 
DGEA4920 
DGEA4930 
DGEA4940 


To ne = - 
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Cerri 


Zo LF 


30 


SS, cONTINUE 


40 


45 


BO 


NN = NO 


((T+HH) .EQ.T) 
write(*,*),’error code = 


KER = 33 
‘ ker 


weep = step + 1 


meate(*,*)’step = 


‘, step 


CALL DGRST (FCN, FCNJ,WK(NY+1) ,WK, WK (NERROR+1) , WK (NSAVE1+4+1) , 


KGO = 


m= 0.DO 


DO 40 I=1,N 


1 WK (NSAVE2+1) ,WK(NPW+1) ,WK (NEQUIL+1) , IWK, NN, step) 


1-KFLAG 
meero (35,55,70,80), 


KGO 
Geen = 1. 23 

NORMAL RETURN FROM INTEGRATOR. THE 
WEIGHTS YMAX(I) ARE UPDATED. IF 
DIFFERENT VALUES ARE DESIRED, THEY 
SHOULD BE SET HERE. A TEST IS MADE 
FOR TOL BEING TOO SMALL FOR THE 
MACHINE PRECISION. ANY OTHER TESTS 
OR CALCULATIONS THAT ARE REQUIRED 
AFTER EVERY STEP SHOULD BE 
INSERTED HERE. IF INDEX = 3, Y IS 
SET TO THE CURRENT SOLUTION ON 
RETURN. IF INDEX = 2, HH IS 
CONTROLLED TO HIT XEND (WITHIN 
ROUNDOFF ERROR), AND THEN THE 
CURRENT SOLUTION IS PUT IN Y ON 
RETURN. FOR ANY OTHER VALUE OF 
INDEX, CONTROL RETURNS TO THE 
INTEGRATOR UNLESS XEND HAS BEEN 
REACHED. THEN INTERPOLATED VALUES 
OF THE SOLUTION ARE COMPUTED AND 
STORED IN Y ON RETURN. 
IF INTERPOLATION IS NOT 
DESIRED, THE CALL TO DGRIN SHOULD 
BE REMOVED AND CONTROL TRANSFERRED 
TO STATEMENT 95 INSTEAD OF 105. 


AYI = DABS (WK (NY+4I) ) 

WK (I) DMAX1 (WK (I) , AYT) 
D = D+(AYI/WK(I))**2 
D = D* (UROUND/TOL) **2 
DN = N 
IF (D.GT.DN) GO TO 75 
IF (INDEX.EQ.3) GO TO 95 
IF (INDEX.EQ.2) GO TO 50 


fo) -XEND) *HH.LT.0.D0) GO TO 25 


NN = NO 


CALL DGRIN 
X = XEND 
oo) 1TO 105 
IF (((T+HH) -XEND) *HH.LE.0O.DO) 


(XEND, WK (NY+1) , NN, Y) 


Sento 25 


IF (DABS (T-XEND) .LE.UROUND*DMAX]1 (10.D0*DABS (T) , HMAX)) GO TO 95 
IF ((T-XEND) *HH.GE.0.D0) GO TO 95 


HH = 
JSTART = 
GO TO 25 


(XEND -T) * (1.D0-4 .DO*UROUND) 


ON AN ERROR RETURN FROM INTEGRATOR, 
AN IMMEDIATE RETURN OCCURS IF 
KFLAG = -2, AND RECOVERY ATTEMPTS 


a7 


DGEA4950 
DGEA49 60 
+ 
DGEA4970 
DGEA49 80 
+ 
et 
DGEA4990 
fe 
DGEA5010 
DGEAS020 
DGEAS 030 
DGEAS040 
DGEAS 050 
DGEAS060 
DGEAS5070 
DGEAS080 
DGEAS090 
DGEA5100 
DGEA5110 
DGEA5120 
DGEA5130 
DGEA5140 
DGEA5150 
DGEA5160 
DGEA5170 
DGEA5180 
DGEA5190 
DGEA5200 
DGEAS5210 
DGEAS5 220 
DGEA5230 
DGEA5240 
DGEA5250 
DGEA5 260 
DGEA5270 
DGEA52 80 
DGEAS 2390 
DGEAS5 300 
DGEA5310 
DGEA5 320 
DGEA5 330 
DGEAS 340 
DGEA5 350 
DGEA5360 
DGEA5 370 
DGEA5 380 
DGEA5 390 
DGEA54 00 
DGEAS410 
DGEAS420 
DGEAS 430 
DGEA5S440 
DGEAS5450 
DGEA5460 
DGEA5470 
DGEAS 4 80 
DGEAS5 490 
DGEAS500 
DGEAS510 
DGEA5520 
DGEAS530 


OAGa 


55 
60 


65 


70 


VS 


80 


85 


96 


75 


100 
205 


110 
9000 


9005 


ARE MADE OTHERWISE. TO RECOVER, HH 
AND HMIN ARE REDUCED BY A FACTOR 
OF .1 UP TO 10 TIMES BEFORE GIVING 


UP. 
JER = 66 
IF (NHCUT.EQ.10) GO TO 65 
MACUL VSaNHeU Tea 
HMIN = HMIN*.1D0 
ft = HE* Spo 
JSTART = -1 
So TO 25 


IF (JER.EQ.66) JER 
IF (JER.EQ.67) JER 
GO TO 95S 


as 2 
136 


Ue = 4 
GO TO 95 


JER = 134 
KFLAG = -2 


so) TO 


JER = 
GO TO 


pe 


67 
60 


WER = 145 
GO TO 1160 


JER = 136 

NN = NO 

CALL DGRIN (XEND, WK (NY+1) , NN, Y) 
X = XEND 

GO TO 110 


= - 

DO 100 Teas aN 

¥ (I) = WK (NY+TI) 

Be (JER 42s) NDEs — KF LAG 
TOUTP = X 

IF (KFLAG.EQ.0) H = HUSED 

IF (KFLAG .NE.0) H = HH 

IER = MAX0 (KER, JER) 
CONTINUE 

IF (KER.NE.O.AND.JER.LT.128) CALL UERTST (KER, 6HDGEAR ) 
tee (JER INE.0) GALT UERTST (JER, 6HDGEAR ) 

RETURN 

END 
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IMSL ROUTINE NAME 


and 


"diag.dat" 


COMPUTER 


LATEST REVISION 


PURPOSE 


PRECISION/HARDWARE 


REQD. IMSL ROUTINES 


NOTATION 


COPYRIGHT 


WARRANTY 


iL 


SUBROUTINE DGRST 


INTEGER 
DOUBLE PRECISION 


= BGEST DGRS0010 
DGRSO0020 
-modified to print sheath and diagnostic output to files "sSheatha.dat" + 
+ 
- IBM/DOUBLE DGRS0050 
DGRSO0060 
= JUNES1, 1982 DGRS0070 
DGRS0080 
- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR DGRS0090 
DGRSO100 
- SINGLE AND DOUBLE/H32 DGRSO0O110 
- SINGLE/H36,H48,H60 DGRS0120 
DGRS0130 
- DGRCS,DGRPS, LUDATF, LUELMF, LEQT1B, UERTST, DGRS0140 
UGETIO DGRSO150 
DGRS0160 
- INFORMATION ON SPECIAL NOTATION AND DGRS0170 
CONVENTIONS IS AVAILABLE IN THE MANUAL DGRSO0180 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP DGRSO0190 
DGRSO0200 
- 1982 BY IMSL, INC. ALL RIGHTS RESERVED. DGRS0210 
DGRSO0220 
- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN DGRS0230 
APPLIED TO THIS CODE. NO OTHER WARRANTY, DGRSO0240 
EXPRESSED OR IMPLIED, IS APPLICABLE. DGRS0250 
DGRS0260 
Siete = em — = = MH Hm mn mm mw mm mw mem em meee ee te eee eee mee ee eee teeter eters DGRS0270 
DGRS 0280 
(FCN, FCNJ, Y, YMAX, ERROR, SAVE1, SAVE2,PW,EQUIL, DGRS0290 
IPIV,NO, step) + 
SPECIFICATIONS FOR ARGUMENTS DGRS0310 
TPIVaG) , NO DGRS0320 
Y(NO,1),YMAX (1) ,ERROR(1) , SAVE1 (1) ,SAVE2 (1), DGRS0330 
PW(1) ,EQUIL(1) ,eprime, eprime (2) + 
SPECIFICATIONS FOR LOCAL VARIABLES DGRS0350 
N,MF,KFLAG, JSTART, NOUSED, NSTEP,NFE,NJE,NSQ, DGRS0360 


PWN PP 


WN 


Oran UP WNP 


INTEGER 


REAL 
DOUBLE PRECISION 


I,METH,MITER,NQ,L, IDOUB, MFOLD, NOLD, IRET, MEO, DGRS0370 
MIO, IWEVAL, MAXDER, LMAX, IREDO,J,NSTEPJ,J1,J2, DGRS0380 
M, IER, NEWO, NPW, NERROR, NSAVE1,NSAVE2,NEQUIL, NY, DGRS0O390 
MITER1, IDUMMY (2) ,NLC,NUC, NWK, JER DGRS0400 

TQ (4) DGRS0410 
T,H, HMIN, HMAX,EPS,UROUND,HUSED,EL(13),OLDLO, DGRS0420 


TOLD, RMAX, RC, CRATE, EPSOLD, HOLD, FN,EDN,E,EUP, DGRS0430 
BND, RH,R1,CON,R,HLO,RO,D, PHLO, PR3,D1,ENQ3, ENQ2, DGRS0440 
PR2,PR1,ENQ1,EPSJ, DUMMY, tcum + 
EXTERNAL FCN, FCNJ DGRS0460 
COMMON /DBAND/ NLC, NUC DGRS0470 
COMMON /GEAR/ T,H,HMIN, HMAX, EPS, UROUND,EPSJ,HUSED, DGRS04 80 
EL, OLDLO, TOLD, RMAX, RC, CRATE,EPSOLD,HOLD, FN, DGRS0490 
EDN,E,EUP,BND,RH,R1,R,HLO,RO,D, PHLO, PR3,D1, DGRSO500 
ENQ3,ENQ2,PR2,PR1,ENQ1,DUMMY,TOQ, DGRS0510 
N,MF,KFLAG, JSTART,NSQO, NOUSED,NSTEP,NFE,NJE, DGRS0520 
NPW, NERROR, NSAVE1, NSAVE2 , NEQUIL, NY, DGRS0530 
I,METH, MITER, NO, L, IDOUB,MFOLD,NOLD, IRET,MEO, DGRS0540 
MIO, IWEVAL, MAXDER, LMAX, IREDO,J,NSTEPJ,J1,dJ2, DGRSO0550 
M, NEWQO, IDUMMY DGRSO0560 
FIRST EXECUTABLE STATEMENT DGRS0570 
open (unit=8,file=’ sheatha.dat’, status=’ unknown’ ) + 
open (unit=9,file=’diag.dat’,status=’ unknown’ ) + 
KFLAG = 0 DGRS0580 
TOD = T DGRS0590 


THIS ROUTINE PERFORMS ONE STEP OF DGRSO0600 


5) 


One) 


Ch. Cl.Cr) ©) 10) 


GG GG) GG) OOOO) O.G@) aman) 
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IF (JSTART GT.0) "GO TO 50 


IF (JSTART.NE.0) 


GO TOL 


THE INTEGRATION OF AN INITIAL 
VALUE PROBLEM FOR A SYSTEM OF 
ORDINARY DIFFERENTIAL EQUATIONS. 


ON THE FIRST CALL, THE ORDER IS SET 


TO 1 AND THE NGTIAL YDOTers 
CALCULATED. RMAX IS THE MAXIMUM 
RATIO BY WHICH H CAN BE INCREASED 
IN A SINGLE STEP? iT IS INITIALLY 
1.E4 TO COMPENSATE FOR THE SMALL 
INITIAL H, BUT THEN IS NORMALLY 
EQUAL TO 10. IF A FAILURE OCCURS 
(IN CORRECTOR CONVERGENCE OR ERROR 
TEST), RMAX IS SET AT 2 FOR THE 
NEXT INCREASE. 


CALL FCN (N,T,Y,SAVE1,eprime, eprime2) 


Bees I=1,N 


Y(I,2) = H*SAVE1 (I) 
METH = MF/10 

MITER = MF-10*METH 
NO = 1 

L= 2 

TDOUB ==. 

RMAX = 1.D4 

RC = 0.D0 

CRATE = 1.D0 

HOLD = “H 

MFOLD = MF 

NSTEP = 0 

NSTEPJ = 0 

NFE = 1 

NJE = 0 

IRET = 3 

GO TO 15 

IF (MF.EQ.MFOLD) GO TO 25 
MEO = METH 

MIO = MITER 

METH = MF/10 

MITER = MF-10*METH 
MFOLD = ME 


IF (MITER.NE.MIO) 


IF (METH.EQ.MEO) 
IDOUB = L+1 
IRET = 1 


IWEVAL 


GO" TO Zs 


IF THE CALLER HAS CHANGED METH, 


DGRCS IS CALLED TO SET THE 
COEFFICIENTS OF THE METHOD. IF THE 
CALLER HAS CHANGED N, EPS, OR 
METH, THE CONSTANTS E, EDN, EUP, 
AND BND MUST BE RESET. EIS A 
COMPARISON FOR ERRORS OF THE 
CURRENT ORDER NO. EUP IS TO TEST 
FOR INCREASING THE ORDER, EDN FOR 
DECREASING THE ORDER. BND IS USED 
TO TEST FOR CONVERGENCE OF THE 
CORRECTOR ITERATES. IF THE CALLER 
HAS CHANGED H, Y MUST BE RESCALED. 
IF H OR METH HAS BEEN CHANGED, 
IDOUB IS RESET TO L + 1 TO PREVENT 
FURTHER CHANGES IN H FOR THAT MANY 
STEPS. 


MITER 


60 


DGRS0610 
DGRS0620 
DGRS0630 


DGRS0640 


DGRSO0650 
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CALL DGRCS (METH,NQ,EL,TQ,MAXDER) 
LMAX = MAXDER+1 


RC = RC*EL(1) /OLDLO 
OLDLO = EL(1) 
FN = N 
EDN = FN*(TQ(1) *EPS) **2 
E = FN*(TQ(2) *EPS) **2 
EUP = FN*(TQ(3) *EPS) **2 
BND = FN* (TQ (4) *EPS) **2 
EPSOLD = EPS 
NOLD = N 
GO TO (30,35,50), IRET 
IF ((EPS.EQ.EPSOLD) .AND. (N.EQ.NOLD)) GO TO 30 
IF (N.EQ.NOLD) IWEVAL = MITER 
ITRET = 1 
GO TO 20 
IF (H.EQ.HOLD) GO TO 50 
RH = H/HOLD 
H = HOLD 
MREDO = 3 
GO TO 40 
RH = DMAX1 (RH, HMIN/DABS (H) ) 
RH = DMIN1 (RH, HMAX/DABS (H) , RMAX) 
i = 1.D0 
DO 45 J=2,L 
Rl = R1*RH 
DO 45 I=1,N 
for, J) = Y(I,J)*R1 
H = H*RH 
RC = RC*RH 
IDOUB = Ltl 


IF (IREDO.EQ.0) GO TO 285 


THiSeoeereon COMPUTES THE PREDICPED 
VALUES BY EFFECTIVELY MULTIPLYING 
THE Y ARRAY BY THE PASCAL TRIANGLE 
MATRIX. RC IS THE RATIO OF NEW TO 


OLD VALUES OF THE COEFFICIENT 


He) 


WHEN RC DIFFERS FROM 1 BY 


MORE THAN 30 PERCENT, OR THE 


CALLER HAS CHANGED MITER, 


TS) SH elOmia@rER TO FORCE THE 
PARTIALS TO BE UPDATED, IF 


PARTIALS ARE USED. 
THE PARTIALS ARE UPDATED AT LEAST 
EVERY 20-TH STEP. 


IF (DABS (RC-1.D0) .GT.0.3D0) IWEVAL = MITER 
IF (NSTEP.GE.NSTEPJ+20) IWEVAL = MITER 
T = T+H 
DO 55 J1=1,NO 
mess O2=J1,NO 
J = (NQ+d1) -J2 
BO 55° I1=1,N 
mete) = Y(I,J)+Y(I,J+1) 


IN ANY CASE, 


UP TO 3 CORRECTOR ITERATIONS ARE 


TAKEN. A CONVERGENCE TEST IS MADE 


ON THE R.M.S. NORM OF EACH 


CORRECTION, USING BND, WHICH IS 
DEPENDENT ON EPS. THE SUM OF THE 
CORRECTIONS IS ACCUMULATED IN THE 


VECTOR ERROR(I). 


NOT ALTERED IN THE CORRECTOR LOOP. 


THE Y ARRAY IS 


ane OPDATED Y VECTOR IS STORED 
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DGRSI230 
DGRS1240 
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DGRS1260 
DGRSY270 
DGRS1280 
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DGRS1300 
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DGRS1440 
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DGRS14 70 
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DGRS1510 
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DGRS1700 
DGRS1710 
DGRS1720 
DGRS1730 
DGRS1740 
DGRS1T 750 
DGRS1760 
DGRS1770 
DGRS1780 
DGRS1790 
DGRS1800 
DGRS1810 
DGRS1820 
DGRS1830 
DGRS1840 
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TEMPORARILY IN SAVE1. DGRS1850 


608be 65 T=i,N DGRS1860 

65 ERROR(I) = 0.DO DGRS1870 

M = 0 DGRS1880 
CALL FCN (N,T,Y,SAVE2,eprime, eprime2) + 

NFE = NFE#+1L DGRS1900 

IF (IWEVAL.LE.0) GO TO 95 DGRS1910 

IF INDICATED, THE MATRIX P = I - DGRS1920 


H*EL(1)*J IS REEVALUATED BEFORE DGRS1930 
STARTING THE CORRECTOR ITERATION. DGRS1940 
IWEVAL IS SET TO 0 AS AN INDICATOR DGRS1950 
THAT THIS HAS BEEN DONE. IF MITER DGRS1960 





= 1OR 2, P IS COMPUTED AND DGRS1970 
PROCESSED IN PSET. IF MITER = 3, DGRS19 80 
THE MATRIX USED IS P =I - DGRS1990 
H*EL(1)*D, WHERE D IS A DIAGONAL DGRS2000 
MATRIX. DGRS2010 
IWEVAL = 0 DGRS2020 
RC = 1.D0 DGRS2030 
NJE = NJE+1 DGRS2040 
NSTEPJ = NSTEP DGRS2050 
GO TO (75,70,80), MITER DGRS2060 
70 NFE = NFE4+N DGRS2070 
75 CON = -H*EL(1) DGRS2080 
MITER1 = MITER DGRS2090 
CALL DGRPS (FCN, FCNJ,Y,N0O,CON,MITER1, YMAX, SAVE1, SAVE2,PW,EQUIL, DGRS2100 
1 IPIV, IER) DGRS2110 
IF (IER.NE.0O) GO TO 155 DGRS2120 
GO TO 125 DGRS2130 
80 R = EL(1)*.1D0 DGRS2140 
BO 85 I=1,N DGRS2150 
8S PW(I) = Y(I,1)+R* (H*SAVE2 (I) -Y(I,2)) DGRS2160 
CALL FCN (N,T, PW, SAVE1,eprime,eprime2) + 
NFE = NFE+1 DGRS2180 
HLO = H*EL(1) DGRS2190 
DO 90 I=1,N DGRS2200 
RO = H*SAVE2 (I) -Y(I, 2) DGRS2210 
PW(I) = 1.D0 DGRS2220 
D = .1D0*RO-H* (SAVE1 (I) -SAVE2 (I) ) DGRS2230 
SAVELGE) = 307 D0 DGRS2240 
IF (DABS (RO) .LT.UROUND*YMAX(I)) GO TO 90 DGRS2250 
IF (DABS (D) .EQ.0.D0) GO TO 155 DGRS2260 
PW(I) = .1D0*RO/D DGRS2270 
SAVE1(I) = PW(I)*RO DGRS2280 
90 CONTINUE DGRS2290 
GO TO 135 DGRS2300 
95 IF (MITER.NE.0) GO TO (125,1257305))) MirTee DGRS2310 | 
DGRS2320 


IN THE CASE OF FUNCTIONAL ITERATION, DGRS2330 
UPDATE Y DIRECTLY FROM THE RESULT DGRS2340 


OF THE LAST FCN CALL. DGRS2350 
Pe= OLD DGRS2360 | 
DO WOO T=1,N DGRS 2370 
R = HeSAVE2 (1) yi 2) DGRS23 80 
D = D+((R-ERROR(I) ) /YMAX(TI) ) **2 DGRS2390 
SAVE1(I) = Y(I,1)+EL(1)*R DGRS2400 
TOG 7ERROR(T) = 7k DGRS2410 
GO TO 145 DGRS 2420 
IN THE CASE OF THE CHORD METHOD, DGRS2430 


COMPUTE THE CORRECTOR ERROR, F SUB DGRS2440 
(M), AND SOLVE THE LINEAR SYSTEM DGRS2450 
WITH THAT AS RIGHT-HAND SIDE AND P DGRS2460 
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AS 


160 
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AS COEFFICIENT MATRIX, USING THE 
LU DECOMPOSITION IF MITER = 1 OR 
2. IF MITER = 3, THE COEFFICIENT 
H*EL(1) IN P IS UPDATED. 
PHLO = HLO 
HLO = H*EL(1) 
IF (HLO.EQ.PHLO) GO TO 115 
R = HLO/PHLO 
DO 110 I=1,N 
D = 1.D0-R*(1.D0-1.D0/PW(I) ) 
IF (DABS(D) .EQ.0.D0) GO TO 165 
PW(I) = 1.D0/D 
DO 120 I=1,N 
SAVE1 (I) = PW(I) * (H*SAVE2 (I) - (Y(I, 2) +ERROR(I)) ) 
GO TO 135 
DO 130 I=1,N 
SAVE1 (I) = H*SAVE2 (I) - (Y(I,2)+ERROR(I) ) 
IF (NLC .EQ. -1) GO TO 131 
NWK = (NLC+NUC+1) *NO41 
CALL LEQT1B(PW,N,NLC,NUC,NO,SAVE1,1,N0,2, PW(NWK) , JER) 
GO TO 135 
CALL LUELMF (PW, SAVE1,IPIV,N,NO, SAVE1) 
D = 0.D0O 
DO 140 I=1,N 
ERROR(I) = ERROR(I)+SAVE1 (I) 
D = D+(SAVE1 (I) /YMAX (I) ) **2 
SAVE1 (I) = Y(1I,1)+EL(1) *ERROR(I) 
TEST FOR CONVERGENCE. IF M.GT.0, THE 
SQUARE OF THE CONVERGENCE RATE 
CONSTANT IS ESTIMATED AS CRATE, 
AND THIS IS USED IN THE TEST. 
IF (M.NE.0) CRATE = DMAX1(.9D0*CRATE,D/D1) 
IF ((D*DMIN1(1.D0,2.D0*CRATE) ) .LE.BND) GO TO 170 
Dl = D 
M = M+l 
IF (M.EQ.3) GO TO 150 
CALL FCN (N,T,SAVE1,SAVE2,eprime,eprime2) 
GO TO 95 
THE CORRECTOR ITERATION FAILED TO 
CONVERGE IN 3 TRIES. IF PARTIALS 
ARE INVOLVED BUT ARE NOT UP TO 
DATE, THEY ARE REEVALUATED FOR THE 
NEXT TRY. OTHERWISE THE Y ARRAY IS 
RETRACTED TO ITS VALUES BEFORE 
PREDICTION, AND H IS REDUCED, IF 
POSSIBLE. IF NOT, A NO-CONVERGENCE 
EXIT IS TAKEN. 
NEE = NFE+2 
IF (IWEVAL.EQ.-1) GO TO 165 
T = TOLD 
RMAX = 2.D0 


mer i60 Jl=1,NO 
moO 160 J2=J1,NO 

J = (NQ+J1) -J2 
mo 160 I=1,N 
meebo) = Y(1,J) -Y(I,J+1) 
IF (DABS (H) .LE.HMIN*1.00001D0) GO TO 280 
Mme= .25D0 
meeDO = 1 
Go TO 35 
IWFVAL = 
GO TO 60 
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THE CORRECTOR HAS CONVERGED. IWEVAL 
iS SET TO -1 IF PARTIAL 
DERIVATIVES WERE USED, TO SIGNAL 
THAT THEY MAY NEED UPDATING ON 
SUBSEQUENT STEPS. THE ERROR TEST 
IS MADE AND CONTROL PASSES TO 
STATEMENT 190 IF IT FAILS. 


IF (MITER.NE.0) IWEVAL = -1l 
NFE = NFE+M 
pa= 0.D0 


DO 175 I=1,N 

D = D+(ERROR(I) /YMAX (I) ) **2 

ie (D.GT.E) GO TOs1g6 

AFTER A SUCCESSFUL STEP, UPDATE THE 

Y ARRAY. CONSIDER CHANGING H IF 
IDOUB = 1. OTHERWISE DECREASE 
IDOUB BY 1. IF IDOUB’'IS THEN 1 AND 
NQ .LT. MAXDER, THEN ERROR IS 
SAVED FOR USE IN A POSSIBLE ORDER 
INCREASE ON THE NEXT STEP. IFA 
CHANGE IN H IS CONSIDERED, AN 
INCREASE OR DECREASE IN ORDER BY 
ONE IS CONSIDERED ALSO. A CHANGE 
IN H IS MADE ONLY IF IT IS BY A 
FACTOR OF AT LEAST 1.1. IF NOT, 


IDOUB IS SET TO 10 TO PREVENT 
TESTING FOR THAT MANY STEPS. 
KFLAG = 0 
IREDO = 0 
NSTEP = NSTEP+1 
PRUSED =) 


NOUSED = NQ 

be 180) J21 a 

DO 180 I=1,N 

Y¥(I,J) = Y¥(I,J)+EL(J) *ERROR (I) 

ir) (LDOUBL HO. 1). Go: Tomo 

IDOUB = IDOUB-1 

IF (IDOUB.GT.1) GO TO 290 

IF (L.EQ.LMAX) GO TO 290 

DO 185 I=1,N 

Y(I,LMAX) = ERROR(I) 

GO TO 290 

THE ERROR TEST FAILED. KFLAG KEEPS 

TRACK OF MULTIPLE FAILURES. 
RESTORE T AND THE Y ARRAY TO THEIR 
PREVIOUS VALUES, AND PREPARE TO 
TRY THE STEP AGAIN. COMPUTE THE 
OPTIMUM STEP SIZE FOR THIS OR ONE 
LOWER ORDER. 

KELAG = KELAG=? 1 

T = TOLD 

DO 195 J1l=1,NO 

DO 195 J2=J1,NQ 

J = (NO today de 

DO 195 I=1,N 

(Te) = YC res ve uol 

RMAX = 2.D0 

IF (DABS (H) .LE.HMIN*1.00001D0) GO TO 270 

ir (KPUAG. EE. -s)8Go TO 260 

MREDO = 2 

Eas = Lope 

GO TO 210 
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DGRS32 80 
DGRS3290 
DGRS3300 
DGRS3310 
DGRS3320 
DGRS3330 
DGRS3340 
DGRS3350 
DGRS3360 
DGRS3370 
DGRS3380 
DGRS3390 
DGRS3400 
DGRS3410 
DGRS3420 
DGRS3430 
DGRS3440 
DGRS3450 
DGRS3460 
DGRS3470 
DGRS3480 
DGRS3490 
DGRS3500 
DGRS3510 
DGRS3520 
DGRS3530 
DGRS3540 
DGRS3550 
DGRS3560 
DGRS3570 
DGRS3580 
DGRS3590 
DGRS3600 
DGRS3610 
DGRS3620 
DGRS3 630 
DGRS3640 
DGRS3650 
DGRS3660 
DGRS3670 
DGRS3680 
DGRS3 690 
DGRS3700 


UXO OYVr QA OA - 


Creare Cr Cree) 


200 


o> 


210 


mes 


220 


22° 


250 


ES > 


240 


245 


2150 


255 


REGARDLESS OF THE SUCCESS OR FAILURE 
OF THEpotlEP, FACTORS PR1, PR2Z, AND 


PR3 ARE COMPUTED, BY WHICH H COULD 


BE DIVIDED AT ORDER NQ - 1, ORDER 


NQ, OR ORDER NO + 1, 


RESPECTIVELY. 


IN THE CASE OF FAILURE, PR3 = 
1.E20 TO AVOID AN ORDER INCREASE. 
THE SMALLEST OF THESE IS 
DETERMINED AND THE NEW ORDER 


CHOSEN ACCORDINGLY. 


TP THE ORDER 


IS TO BE INCREASED, WE COMPUTE ONE 
ADDITIONAL SCALED DERIVATIVE. 


meso = 1.D+20 
me (L.EQO.LMAX) GO TO 210 
D1 = 0.DO 


DO 205 I=1,N 

D1 = D1+((ERROR(I) -Y (I, LMAX) ) /YMAX(I)) **2 
ENQ3 = .5D0/(L+1) 

PR3 = ((D1/EUP) **ENQ3) *1.4D0+1.4D-6 

ENQ2 = .5D0/L 

PR2 ((D/E) **ENQ2) *1.2D0+1.2D-6 


PR1 = 1.D+20 

IF (NQ.EQ.1) GO TO 220 

D = 0.DO 

DO 215 I=1,N 

D = D+(Y¥(I,L) /YMAX (I) ) **2 

ENQ1 = .5D0/NQ 

PR1 = ((D/EDN) **ENQ1) *1.3D0+1.3D-6 
IF (PR2.LE.PR3) GO TO 225 

IF (PR3.LT.PR1) GO TO 235 

GO TO 230 

IF (PR2.GT.PR1) GO TO 230 

NEWQ = NOQ 

RH = 1.D0/PR2 

GO TO 250 

NEWQ = NQ-1 

RH = 1.D0/PR1 

IF (KFLAG.NE.0.AND.RH.GT.1.D0) RH = 1.D0 
GO TO 250 

NEWQ = L 

RH = 1.D0/PR3 

Mee ({RH.LT.1.1D0) GO TO 245 

DO 240 I=1,N 

Y(I,NEWQ+1) = ERROR(I) *EL(L) /L 
GO TO 255 

IDOUB = 10 

GO TO 290 


IF ((KFLAG.EQ.0) .AND. (RH.LT.1.1D0)) GO TO 245 


IF (NEWQ.EQ.NQ) GO TO 35 
NO = NEWO 


L = NQ+1 
Peet = 2 
ee TO 15 


IF THERE IS A CHANGE OF ORDER, RESET 
NO, L, AND THE COEFFICIENTS. IN 
ANY CASE H IS RESET ACCORDING TO 
RH AND THE Y ARRAY IS RESCALED. 
THEN EXIT FROM 285 IF THE STEP WAS 
OK, OR REDO THE STEP OTHERWISE. 


CONTROL REACHES THIS SECTION IF 3 OR 
MORE FAILURES HAVE OCCURED. IT IS 
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DGRS3710 
DGRS3720 
DGRS3730 
DGRS3740 
DGRS3750 
DGRS3760 
DGRS3770 
DGRS3780 
DGRS3790 
DGRS3 800 
DGRS3 810 
DGRS3 820 
DGRS3830 
DGRS3 840 
DGRS3850 
DGRS3 860 
DGRS3870 
DGRS3880 
DGRS3 890 
DGRS3900 
DGRS3910 
DGRS3920 
DGRSs7 50 
DGRS3940 
DGRS3 930 
DGRS3960 
DGRS3970 
DGRS39 80 
DGRS3990 
DGRS4000 
DGRS4010 
DGRS4020 
DGRS4030 
DGRS4040 
DGRS4050 
DGRS4060 
DGRS4070 
DGRS4080 
DGRS4090 
DGRS4100 
DGRS4110 
DGRS4120 
DGRS4130 
DGRS4140 
DGRS4150 
DGRS4160 
DGRS4170 
DGRS4180 
DGRS4190 
DGRS4200 
DGRS4210 
DGRS4220 
DGRS4230 
DGRS4240 
DGRS4250 
DGRS4260 
DGRS4270 
DGRS4280 
DGRS4290 
DGRS4300 
DGRS4310 
DGRS4 320 


COO OOOO GO) 


260 ee 


Z65 


CAO 


Z27e 


275 


280 


285 
290 


RH . 1D 


ee H* RE 


(KFLAG.EO.-7) GO 8O.275 


0 


RH = DMAX1 (HMIN/DABS (H) , RH) 


ASSUMED THAT THE DERIVATIVES THAT 
HAVE ACCUMULATED IN THE Y ARRAY 
HAVE ERRORS OF THE WRONG ORDER. 
HENCE THE FIRST DERIVATIVE IS 
RECOMPUTED, AND THE ORDER IS SET 
TO 1. THEN H IS REDUCED BY A 
FACTOR OF 10, AND THE STEP IS 
RETRIED. AFTER A TOTAL OF 7 
FAILURES, AN EXIT IS TAKEN WITH 
KFEAG = -2. 


CALL FCN (N,T,Y,SAVE1,eprime, eprime2) 


NFE = NFE+1 
HOs265 1=-1,N 
H*SAVE1 (I) 
MITER 

10 
IF (NQ.EQ.1) 


Mir,2) = 
IWEVAL = 
IDOUB = 


NO-= 1 

L = 2 
IRET = 3 
GO TO 15 


NQ 


ALL RETURNS ARE MADE THROUGH THIS 
SECTION. H IS SAVED IN HOLD TO 
ALLOW THE CALLER TO CHANGE H ON 
THE NEAT STEP 


C--Diagnostic Check of first and second derivatives of E 


300 


305 


ac (tcum.eq.told) go to 310 


write (8,300) tcum,step,y(1,1),Vvi2,1) ,v(3,.)) (v4) Vv oe 
format (1x,e11.4,1x,15,5(1x,e11.4)) 
write (9,305) step, eprime,eprime2 


format (1x,1I5,2 (1x,e20.13) ) 


RETURN 
END 
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DGRS4330 
DGRS4 340 
DGRS4350 
DGRS4 360 
DGRS4370 
DGRS4380 
DGRS4390 
DGRS4400 
DGRS4410 
DGRS4420 
DGRS4430 
DGRS4440 
DGRS4450 
DGRS4460 
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DGRS4480 
DGRS4490 
DGRS4500 
DGRS4510 
DGRS4520 
DGRS4530 
DGRS4540 
DGRS4550 
DGRS4560 
DGRS4570 
DGRS4580 
DGRS4590 
DGRS4600 
DGRS4610 
DGRS4620 
DGRS4630 
DGRS4640 
DGRS4650 
DGRS4660 
DGRS4670 
DGRS4680 
DGRS4690 
DGRS4700 

we 
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DGRS4710 
DGRS4720 
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DGRC0010 
DGRC0020 


DGRC0040 
DGRC0050 
DGRC0060 


IMSL ROUTINE NAME - Uekes 

wee ee ne ne ee ee ee ee ee ee ee ee eee DGRC0030 
COMPUTER - IBM/DOUBLE 
LATEST REVISION - JANUARY 1, 1978 


PURPOSE 


PRECISION/HARDWARE 


- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 


- SINGLE AND DOUBLE/H32 


- SINGLE/H36,H48,H60 


REQD. IMSL ROUTINES 


NOTATION - 


- NONE REQUIRED 


INFORMATION ON SPECIAL NOTATION AND 


CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 


COPYRIGHT - 


WARRANTY > 


IMSL WARRANTS ONLY THAT 
MPPLIED TO THIS CODE. 
EXPRESSED OR IMPLIED, 


1978 BY IMSL, INC. ALL RIGHTS RESERVED. 


IMSL TESTING HAS BEEN 
NO OTHER WARRANTY, 
IS APPLICABLE. 


DGRC0070 
DGRC0080 
DGRC0090 
DGRC0100 
DGRC0110 
DGRC0120 
DGRC0130 
DGRC0140 
DGRCOPSO 
DGRC0160 
DGRC0170 
DGRC0180 
DERCOLI0 
DGRC0200 
DGRC0210 
DGRC0220 
DGRC0230 
DGRC0240 
DGRC0250 


ee ee ici =~ ~ = ~ - - - 5 2 ein a ee ee ee ~~~ - +e - DGRC0260 


SUBROUTINE DGRCS 


plas 
-clQ2Z2E-2, -2945n-35, 5492-4, 
COOH = 6) Ae Si 66am 
eel G7b-1, 7°15. ,2.)L2qee os eo oo. 

See on OfO0o,.67.97,10Gee, 126.7, 
Magma, wooo, tol .0,2 7074550 435 34 
MOna2,135.7,7°1.,12.072450 eee, 
Somos OOo no).91;, LOGes,.2o.7 | 
tageray hoo. Oo, 191.0,12 3076.0) 


(METH, NO, EL, TQ, MAXDER) 


SPECIFICATIONS FOR ARGUMENTS 


SPECIFICATIONS FOR LOCAL VARIABLES 


eee oO. — Al. 


INTEGER METH, NO, MAXDER 
REAL TQ (1) 

DOUBLE PRECISION EL (1) 

INTEGER K 

REAL PoRroE (12,2, 3) 
DATA PERTST/1. 
- olla al) oy 
ze "20g 25-9,, 
3 
4 
5 

6 

7 

8 
9 


BeerO (5,10), 
S MAXDER = 12 


METH 


Behera. 57 8~ 1. / 


FIRST EXECUTABLE STATEMENT 


meee, (15,20,25,30,35,40,45,50,55,60,65,70), NO 


10 MAXDER = 5 


meee (75,80,85,90,95), NO 


THE FOLLOWING COEFFICIENTS SHOULD BE 
DEFINED TO MACHINE ACCURACY. FOR A 
GIVEN ORDER NO, THEY CAN BE 
CALCULATED BY USE OF THE 
GENERATING POLYNOMIAL L(T), 
COEFFICIENTS ARE EL(I).. 
EL (1) + ELQ)4Yr =)... + 
EL(NQ+1) *T**NQ. FOR THE IMPLICIT 
ADAMS METHODS, L(T) IS GIVEN BY 
DL/DT = (T+1) * (T+2) * 
* (T+NQO-1)/K, L(-1) = 


WHOSE 
L(T) = 


O, WHERE K = 
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DGRC0270 
DGRC0280 
DGRC0290 
DGRC0300 
DGRC0310 
DGRC0320 
DGRC0330 
DGRC0340 
DGRCO0350 
DGRC0360 
DGRC0370 
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DGRC0600 
DGRC0610 
DGRC0620 
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20 
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Do 


EL (10) 


al 


0 


© 


be 
0 
OT 
OG: 
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om 
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4. 
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© 
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Or. 
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0 
3). 
ie 
Te 
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4. 
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a 
0 
3s 
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ie 
ee 
Le. 
ye 
Zn 
0 
ae 
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4. 
Dee: 
On 
Le. 
i 
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nnuuunwuntwanwreruad but uu i 
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ODO 


5D0 
5D0 


166666666666667D-01 
75D0 
666666666666667D-01 


375D0 

166666666666667D-01 
2333355552355 555-08 
166666666666667D-02 


486111111111111D-01 
041666666666667D0 

B861LLITZI1TII1 111 eo 
041666666666667D-01 
333333333333333De0— 


29S6TTTIT TTT t 
141666666666667D+00 
625D+00 
71708333333333e83Rr0m 
025D+00 
388888888888889D- 03 


1559193221693 12830). 
225D+00 

2LB8518518SresioEpeoL 
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S6lLILTIITT Aaa aos 
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685165185185 185D-02 
357638888888889D-01 
TI777777 tele 7 78D 
064814814814815D-02 
936507936507937D-04 
480158730158730D-05 


948680004409171D-01 
358928571428571D+00 
765542328042328D-01 
1718750-07 
113541666666667D-01 
01875D+00 
934523809523810D-03 
116071428571429D-04 
gan /31IZZ3985389D- Ce 


FACTORIAL (NQ-1) . 
METHODS, L(T) = 
*(T+NQ) /K, WHERE K = 

FACTORIAL (NQ)*(1 + 1/2 + 
1/NQ) . 
GROUPS APPEAR BELOW IS.. 


FOR THE 


THE ORDER IN WHIGH THE 


GEAR 


(T+1) * (T+2) * 


+ 


IMPLICIT 


ADAMS METHODS OF ORDERS 1 TO 12, 
BACKWARD DIFFERENTIATION METHODS 


OF SORDERS) 1 -TOe5: 
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DGRC0630 
DGRC0640 
DGRC0650 
DGRC0660 
DGRC0670 
DGRC0680 
DGRC0690 
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80 
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25 
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GO TO 
EL (1) 
EL (3) 
EL (4) 
EL (5) 
EL (6) 
EL (7) 
EL(8) 
EL (9) 


Bat 20) = 
Pili) <= 


GO TO 
EL (1) 
EL (3) 
EL (4) 
EL (5) 
EL (6) 
EL (7) 
EL (8) 
EL (9) 


(10) = 
BG) = 
EL (12) 


GO TO 
EL (1) 
EL (3) 
EL (4) 
EL (5) 
EL (6) 
EL (7) 
EL (8) 
EL (9) 


EL(10) 
EL (11) 
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.906057098765432D-02 
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.171514550264550D+00 
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-591565255 /omapi224D-06 
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tt th Ua We 


e il 
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-742655400315991D-01 
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.260271164021164D+00 
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ogveto20523222D-02 
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=o Og 1 5 Soe mes44p-05 
= 4-9322530864197531D-06 


tet t wou tf uw ow 


FPUWUNHDRFPHEFPNONKFAEHHAPKFPUPRPNONKFNWNER PHP EN 


mime 2)= 1.503126503126503D-07 
EL(13)= 2.087675698786810D-09 


so TO 


EL (1) 
GO TO 
Be (1) 
EL (3) 
GO TO 
EL (1) 
EL (3) 
EL (4) 
GO TO 
EL (1) 
EL (3) 
EL (4) 
EL(5) 
GO TO 
ea (1) 
EL (3) 
EL (4) 
EL (5) 
EL (6) 


.0D+00 


2 
0 
6.666666666666667D-01 
See 3 3333333338D-01 
0 
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.454545454545455D-01 
EL (1) 
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.48D+00 
. 7D+00 
- 2D+00 
.02D+00 


mee g5o2043795620D-01 
yaretle78832116788D-01 
PeOzreo751021898D-01 
-474452554744526D-02 
-649635036496350D-03 


voutunuielbrePnrunt tf end bed tt eo 
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Boer10S K=1,3 
TOQO(K) = PERTST(NQ, METH, K) 
CONTINUE 


TO (4) 


= .5D0*TO(2) / (NQ+2) 


DGRC1250 
DGRC1260 
DGRC1270 
DGRC1280 
DGRC1290 
DGRC1300 
DGRC1310 
DGRC1320 
PGE 0 
DGRC1340 
DGRC1350 
DGRC1360 
DGRC1370 
DGRC13 80 
DERG 70 
DGRC1400 
DGRC1410 
DGRC1420 
DGRC1430 
DGRC1440 
DGRC1450 
DGRC1460 
DGRC1470 
DGRC1480 
DGRC1490 
DGRC1500 
DGRC1510 
DGRC1520 
DGRC1530 
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RETURN = 
END 
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IMSL ROUTINE NAME -“DGRES DGRPO0O10 


DGRP0020 
eee +--+ ---- DGRP0030 
DGRP0040 

COMPUTER - IBM/DOUBLE DGRPOOSO 
DGRP0060 

LATEST REVISION - NOVEMBER 1, 1984 DGRP0070 
DGRP0O080 

PURPOSE - NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR DGRP0O090 
DGRP0100 

PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 DGRP0O110 
- SINGLE/H36,H48,H60 DGRP0120 

DGRP0130 

REOD. IMSL ROUTINES - LUDATF, LEQT1B, UERTST, UGETIO DGRP0140 
DGRP0O150 

NOTATION - INFORMATION ON SPECIAL NOTATION AND DGRP0160 
CONVENTIONS IS AVAILABLE IN THE MANUAL DGRP0170 

INTRODUCTION OR THROUGH IMSL ROUTINE UHELP DGRP0180 

DGRP0190 

COPYRIGHT - 1984 BY IMSL, INC. ALL RIGHTS RESERVED. DGRP0200 
DGRP0210 

WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN DGRP0220 
APPLIED TO THIS CODE. NO OTHER WARRANTY, DGRP0230 

EXPRESSED OR IMPLIED, IS APPLICABLE. DGRP0240 

DGRP0250 

eae a ee a eee ee DGRP0260 
DGRP0270 

SUBROUTINE DGRPS (FCN, FCNJ,Y,NO,CON,MITER, YMAX,SAVE1,SAVE2,PW, DGRP0O280 

* EQUIL, IPIV, IER) DGRP0290 
SPECIFICATIONS FOR ARGUMENTS DGRP0300 

INTEGER NO,MITER, IPIV(1),IER DGRP0310 
DOUBLE PRECISION Y(NO,1),CON, YMAX(1),SAVE1(1),SAVE2(1),PW(1), DGRP0320 

* EQUIL (1) DGRP0330 
SPECIFICATIONS FOR LOCAL VARIABLES DGRP0340 

DGRP0350 

INTEGER NC,MFC,KFLAG, JSTART, NOUSED,NSTEP,NFE,NJE,NPW, DGRP0360 

* NSQ,1I,J1,70,NERROR, NSAVE1, NSAVE2, NEQUIL, NY, DGRP0370 
* IDUMMY (23) ,NLIM, II, IJ, LIM1, LIM2,NB, NLC, NUC, NWK DGRPO380 
REAL SDUMMY (4) DGRP0390 
DOUBLE PRECISION T,H,HMIN, HMAX,EPSC,UROUND, EPSJ,HUSED,D,RO,YJ,R,DGRP0400 

* D1,D2, WA, DUMMY (40) DGRP0410 
COMMON /DBAND/ NLC, NUC DGRP0420 
COMMON /GEAR/ T,H, HMIN, HMAX, EPSC, UROUND, EPSJ, HUSED, DUMMY, DGRP0430 

* SDUMMY ,NC,MFC,KFLAG, JSTART,NSQ,NQUSED,NSTEP, DGRP0440 
* NFE,NJE,NPW, NERROR, NSAVE1, NSAVE2, NEQUIL, NY, DGRP0450 
* I DUMMY DGRP0460 


THIS ROUTINE IS CALLED BY DGRST TO DGRP0470 
COMPUTE AND PROCESS THE MATRIX P = DGRP0480 
Peon eiei!) *J , WHERE J iS AN DGRP0490 
APPROXIMATION TO THE JACOBIAN. J DGRPO0500 
LS COMBUIED, EDITHER BY THE USER= DGRP0510 
SUPPLIED ROUTINE FCNJ IF MITER = DGRP0520 
1, OR BY FINITE DIFFERENCING IF DGRP0530 
MITER = 2. J IS STORED IN PW AND DGRP0540 
REPLACED BY P, USING CON = DGRPO550 
“Heeb THEN P IS SUBJECTED TO DGRP0560 
LU DECOMPOSITION IN PREPARATION DGRP0570 
FOR LATER SOLUTION OF LINEAR DGRP0580 
SYSTEMS WITH P AS COEFFICIENT DGRP0590 
MATRIX. IN ADDITION TO VARIABLES DGRP0600 
DESCRIBED PREVIOUSLY, DGRP0610 


COMMUNICATION WITH DGRPS USES THE DGRP0620 


ye. 


(11) 


OQ 


10 


5 


Zo 


LD 


30 


35 


40 


45 


50 


5 


60 


FOLLOWING EPSJ = DSQRT(UROUND) , 
USED IN THE NUMERICAL JACOBIAN 
INCREMENTS. 


FIRST EXECUTABLE STATEMENT 

IF (NLC.EQ.-1) GO TO 45 

BANDED JACOBIAN CASE 
NB = NLC+NUC+1 
NWK = NB*NO+1 
IF (MITER.EQ.2) GO TO 15 

MITER = 1 
NLIM = NB*NO 
DO S I=1,NLIM 

PW(I) = 0.0D0 
CONTINUE 
CALL FCNJ (NC,T,Y, PW) 
DO 10 I=1,NLIM 
PW(I) = PW(I) *CON 
CONTINUE 
GO TO 35 
MITER = 2 

D = 0.0D0 
DO® 20 I=1,Ne 
D = D+SAVE2 (I) **2 
RO = DABS (H) *DSQORT (D) *1.0D+03 *UROUND 
pe 30 J=1),NC 


YJ = YC 1 

R = EPSJ*YMAX (J) 
R = DMAX1(R,RO) 
Y(0,1)) =44 Gade 
D = CON/R 


CALL ECN (NG iy SAVE) 
LIM1 = MAXO(1,J-NUC) 
LIM2 = MINO(NO,J+NLC) 
DO 25. Labi Mier 


IJ = (J-I+NLC) *NO+I 
PW(IJ) = (SAVE1 (I) -SAVE2 (I) )*D 
CONTINUE 
Y(J,1) = Yd 
CONTINUE 


ADD IDENTITY MATRIX. 
BO 40 IT=1,NC 
Il =,NLC*Noet 
PW(II) = PW(II)+1.0D0 
CONTINUE 


DO LU DECOMPOSITION ON P 


CALL LEQT1B (PW,NC,NLC,NUC,NO,EQUIL,1,N0,1,PW(NWK) , IER) 
RETURN 


FULL JACOBIAN CASE 
TF (MITERVEO?2Z ) BGO 1e 55 


MITER = 1 
CALL FCNJ (NC,T,Y, PW) 
DO 50. 1.=1.,,Ns0 
PW(I) = PW(I) *CON 
GO TO 75 

MITER = 2 
D = 0.0D0 


BO 60 f= 

D = D+SAVE2 (I) **2 

RO DABS (H) *DSQRT (D) *1 .0D+03*UROUND 
J1 0 


We 


DGRP0630 
DGRP0640 
DGRP0650 
DGRP0660 
DGRP0670 
DGRP0680 
DGRP0690 
DGRP0700 
DGRP0710 
DGRP0720 
DGRP0730 
DGRP0740 
DGRP0750 
DGRP0760 
DGRP0770 
DGRP0780 
DGRP0790 
DGRPO0 800 
DGRP0810 
DGRP0820 
DGRPO 830 
DGRP0840 
DGRPO0850 
DGRPO860 
DGRP0870 
DGRP0880 
DGRP0890 
DGRP0900 
DGRP0910 
DGRP0920 
DGRP0930 
DGRP0940 
DGRP0950 
DGRP0960 
DGRP0970 
DGRPO980 
DGRP0990 
DGRP1000 
DGRP1010 
DGRP1020 
DGRP1030 
DGRP1040 
DGRP1050 
DGRP1060 
DGRP1070 
DGRP1080 
DGRP1090 
DGRP1100 
DGRP1110 
DGRP1120 
DGRP1130 
DGRP1140 
DGRP1150 
DGRP1160 
DGRP1170 
DGRP1180 
DGRP1190 
DGRP1200 
DGRP1210 
DGRP1220 
DGRP1230 
DGRP1240 
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70 


75 


80 


mo 70 J=1,NC 
eel = Ys, a} 
R EPSJ* YMAX (J) 
R = DMAX1 (R, RO) 
mo, 1) = Y(J,1)+R 
D = CON/R 
CALL FCN(NC,T,Y,SAVE1) 
HO 65y1=1,NC 


PW(I+J1) = (SAVE1 (I) -SAVE2(I))*D 
mel) = YJ 
J1l = J1+N0 
CONTINUE 
ADD IDENTITY MATRIX. 
a= 1 


DO 80 I=1,NC 
PW(J) = PW(J)+1.0D0 
J = J+(NO+1) 
CONTINUE 


DO LU DECOMPOSITION ON P. 


CALL LUDATF (PW, PW,NC,NO,0,D1,D2, IPIV, EQUIL, WA, IER) 
RETURN 
END 


73 


DGRP1250 
DGRP1260 
DGRP1270 
DGRP1280 
DGRP1290 
DGRP1300 
DGRP1310 
DGRP1320 
DGRP1330 
DGRP1340 
DGRP1350 
DGRP1360 
DGRP1 370 
DGRP13 80 
DGRP1390 
DGRP1400 
DGRP1410 
DGRP1420 
DGRP1430 
DGRP1440 
DGRP1450 
DGRP1460 
DGRP1470 
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IMSL ROUTINE NAME 


DGRIO010 
DGRI0020 


ween ee ee ee re re re ee ee ee ee ee ee +--+ --- DGRI0030 


COMPUTER 
LATEST REVISION 
PURPOSE 


PRECISION/HARDWARE 


REQD. IMSL ROUTINES 


NOTATION 


GOPYRIGHT 


WARRANTY 


- IBM/DOUBLE 


- JANUARY 1, 1978 


- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 


- SINGLE AND DOUBLE/H32 
- SINGLE/H36,H48,H60 


- NONE REQUIRED 


- INFORMATION ON SPECIAL NOTATION AND 
CONVENTIONS IS AVAILABLE IN THE MANUAL 


INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 


- 1978 BY IMSL, INC. ALL RIGHTS RESERVED. 


- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 


APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


DGRI0040 
DGRIOOS0O 
DGRIO0060 
DGRI0070 
DGRI0080 
DGRIO090 
DGRI0100 
DGRI0110 
DGRI0120 
DGRI0130 
DGRI0140 
DERIUCiSo0 
DGRI0160 
DGRI0170 
DGRI0180 
DGRI0190 
DGRIO0200 
DGRI0210 
DGRI0220 
DGRI0230 
DGRI0240 
DGRIO0250 


ee en en ee ee ee ee ee ee ee ee eee DGRI0260 


SUBROUTINE DGRIN 


INTEGER 
DOUBLE PRECISION 


INTEGER 
iL 
2 

REAL 

DOUBLE PRECISION 
al 

COMMON /GEAR/ 


i 
2 
3 


DOS alee NG 
OnE) = ¥ Cl) 
5 CONTINUE 


JSTART + 1 

(TOUT - T)/H 

S1 = 1.0D0 

DO i15s0)——2, 
‘Sl = S1*S 


L 
S 


(TOUT, Y, NO, YO) 

SPECIFICATIONS FOR ARGUMENTS 
NO 
TOUT, YO (NO) , Y(NO,1) 

SPECIFICATIONS FOR LOCAL VARIABLES 
NC,MFC,KFLAG,1I,L,J,JSTART, NSQ, NQUSED, NSTEP, 
NFE,NJE, NPW, NERROR, NSAVE1,NSAVE2,NEQUIL,NY, 

IDUMMY (23) 

SDUMMY (4) 

T,H, HMIN, HMAX, EPSC, UROUND, EPSJ, HUSED,S,S1, 
DUMMY (40) 

T,H,HMIN, HMAX, EPSC, UROUND, EPSJ, HUSED, DUMMY, 
SDUMMY , NC, MFC, KFLAG, JSTART, NSQ, NQUSED,NSTEP, 
NFE, NUE, NPW, NERROR, NSAVE1, NSAVE2, NEQUIL,NY, 
I DUMMY 

FIRST EXECUTABLE STATEMENT 


DGRI0270 
DGRI0280 
DGRI0290 
DGRI0300 
DGRI0310 
DGRI0320 
DGRI0330 
DGRI0340 
DGRIO0350 
DGRI0360 
DGRI0370 
DGRI0380 
DGRI0390 
DGRI0400 
DGRI0410 
DGRI0420 
DGRI0430 
DGRI0440 
DGRI0450 
DGRI0460 


THIS SUBROUTINE COMPUTES INTERPOLATEDDGRI0470 


VALUES OF THE DEPENDENT VARIABLE 
Y AND STORES THEM IN YO. THE 
INTERPOLATION IS TO THE 

POINT TY =" TOUL, AND USES, THE 
NORDSIECK HISTORY ARRAY Y, AS 


FOLLOWS. . 
NQ 

YO (2) = SUM Y(I,J+1) *S**J 
J=0 

WHERE S = - (T-TOUT) /H. 
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DGRI0480 
DGRI0490 
DGRIO500 
DGRI0510 
DGRI0520 
DGRI0530 
DGRI0540 
DGRIO550 
DGRIO0560 
DGRI0570 
DGRIO580 
DGRI0590 
DGRI0600 
DGRI0610 
DGRI0620 


<< 


——— oo 





pesto. = 1,NC 


YO (TI) 
CONTINUE 


==YO(l) + <Siey(t, 7) 


eo 


DGRI0630 
DGRI0640 
DGRI0650 
DGRI0660 
DGRI0670 
DGRI0680 


0A 


AN MAA ANDADONAANNNAANNANNDANANAANNADAAARAAADAARDNDAAADAADN1AAH1O21An Oe ee a 


IMSL ROUTINE NAME 


COMPUTER 


LATEST REVISION 


PURPOSE 


USAGE 


ARGUMENTS 


A 


LU 


LA 


IDGT 


me 


D2 


IPVT 


EQUIL 


WA 


IER 


- LUDATF 


LUDA0010 
LUDA0020 
LUDA0030 


LUDAO0040 


- IBM/DOUBLE 


- JANUARY 1, 1978 
- L-U DECOMPOSITION BY THE CROUT ALGORITHM 
WITH OPTIONAL ACCURACY TEST. 


- CALL) LUPETE (Ay nUnN, In, LDGT,D1,D2,5Pun 
EQUIL, WA, IER) 


- INPUT MATRIX OF DIMENSION N BY N CONTAINING 
THE MATRIX TO BE DECOMPOSED. 

- REAL OUTPUT MATRIX OF DIMENSION N BY N 
CONTAINING THE L-U DECOMPOSITION OF A 
ROWWISE PERMUTATION OF THE INPUT MATRIX. 
FOR A DESCRIPTION OF THE FORMAT OF LU, SEE 
EXAMPLE. 

- INPUT SCALAR CONTAINING THE ORDER OF THE 
MATRIX A. 

- INPUT SCALAR CONTAINING THE ROW DIMENSION OF 
MATRICES A AND LU EXACTLY AS SPECIFIED IN 
THE CALLING PROGRAM. 

INPUT OPTION. 

IF IDGT IS GREATER THAN ZERO, THE NON-ZERO 
ELEMENTS OF A ARE ASSUMED TO BE CORRECT TO 
IDGT DECIMAL PLACES. LUDATF PERFORMS AN 
ACCURACY TEST TO DETERMINE IF THE COMPUTED 
DECOMPOSITION IS THE EXACT DECOMPOSITION 
OF A MATRIX WHICH DIFFERS FROM THE GIVEN 
ONE BY LESS THAN ITS UNCERTAINTY. 

IF IDGT IS EQUAL TO ZERO, THE ACCURACY TEST 
IS BYPASSED. 

- OUTPUT SCALAR CONTAINING ONE OF THE TWO 
COMPONENTS OF THE DETERMINANT. SEE 
DESCRIPTION OF PARAMETER D2, BELOW. 

- OUTPUT SCALAR CONTAINING ONE OF THE 
TWO COMPONENTS OF THE DETERMINANT. THE 
DETERMINANT MAY BE EVALUATED AS (D1) (2**D2). 

- OUTPUT VECTOR OF LENGTH N CONTAINING THE 
PERMUTATION INDICES. SEE DOCUMENT 
(ALGORITHM) . 

- OUTPUT VECTOR OF LENGTH N CONTAINING 
RECIPROCALS OF THE ABSOLUTE VALUES OF 
THE LARGEST (IN ABSOLUTE VALUE) ELEMENT 
IN EACH ROW. 

- ACCURACY TEST PARAMETER, OUTPUT ONLY IF 
IDGT IS GREATER THAN ZERO. 

SEE ELEMENT DOCUMENTATION FOR DETAILS. 
- ERROR PARAMETER. (OUTPUT) 
TERMINAL ERROR 
IER = 129 INDICATES THAT MATRIX A IS 
ALGORITHMICALLY SINGULAR. (SEE THE 
CHAPTER L PRELUDE) . 
WARNING ERROR 
IER = 34 INDICATES THAT THE ACCURACY TEST 
FAILED. THE COMPUTED SOLUTION MAY BE IN 
ERROR BY MORE THAN CAN BE ACCOUNTED FOR 
BY THE UNCERTAINTY OF THE DATA. THIS 
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LUDAO0050 
LUDA0060 
LUDA0070 
LUDA0080 
LUDA0090 
LUDA0100 
LUDA0110 
LUDA0N120 
LUDA0130 
LUDA0140 
LUDA0150 
LUDA0160 
LUDA0170 
LUDA0180 
LUDA0190 
LUDA0200 
LUDA0210 
LUDA0220 
LUDA0230 
LUDA0240 
LUDAO0250 
LUDA0260 
LUDA0270 
LUDAO02 80 
LUDA0290 
LUDA03 00 
LUDA0310 
LUDA0320 
LUDA0330 
LUDA0340 
LUDA0350 
LUDA03 60 
LUDA0370 
LUDA03 80 
LUDA039 0 
LUDA0400 
LUDA0410 
LUDA0420 
LUDA04 30 
LUDA0440 
LUDA0450 
LUDA0460 
LUDA0470 
LUDA04 80 
LUDA0490 
LUDAO0500 
LUDA0510 
LUDA0520 
LUDA0530 
LUDA0540 
LUDA0550 
LUDA0560 
LUDAO0S5 70 
LUDAOS5 80 
LUDAO0590 
LUDA0600 
LUDA0610 
LUDA0 620 


OO OY COA ae eae eee 0.) 


o) 


PRECISION/HARDWARE 


REQD. 


NOTATION 


REMARKS 


eer y RIGHT - 


WARRANTY 


IMSL ROUTINES 


WARNING CAN BE PRODUCED ONLY IF IDGT IS 
GREATER THAN 0 ON INPUT. SEE CHAPTER L 
PRELUDE FOR FURTHER DISCUSSION. 


- SINGLE AND DOUBLE/H32 
- SINGLE/H36,H48,H60 


= UERTST, UGETIO 


- INFORMATION ON SPECIAL NOTATION AND 
CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 


A TEST FOR SINGULARITY IS MADE AT TWO LEVELS: 
1. A ROW OF THE ORIGINAL MATRIX A IS NULL. 
2. A COLUMN BECOMES NULL IN THE FACTORIZATION PROCESS. 


1378 BY IMSL, INC. ALL RIGHTSsRESERVED: 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


SUBROUTINE LUDATF (A,LU,N,IA, IDGT,D1,D2,IPVT, EQUIL, WA, IER) 


DIMENSION 


Paha 1), LU Pa eben (1) EOULL (1) 


DOUBLE PRECISION A,LU,D1,D2,EQUIL, WA, ZERO, ONE, FOUR, SIXTN, SIXTH, 


* 


RN, WREL, BIGA, BIG, P, SUM,AI,WI,T, TEST,Q 


DATA ZERO, ONE, FOUR, SIXTN, SIXTH/0.D0,1.D0,4.D0, 
* 16.D0>. 0625D07 
FIRST EXECUTABLE STATEMENT 
INITIALIZATION 
IER = 0 
RN = 
WREL = ZERO 
Dl = ONE 
D2 = ZERO 
BIGA = ZERO 
moO 10 I=1,N 
BIG = ZERO 
BO S J=1,N 
P = A(I,J) 
EU (1, 0). = 
P = DABS (P) 
IF (P .GT. BIG) BIG = P 
5 CONTINUE 
IF (BIG .GT. BIGA) BIGA = BIG 
mE (BIG .BO. ZERO) GO TO 110 
EQUIL(I) = ONE/BIG 
10 CONTINUE 
moO 105 J=1,N 
wml = J-1 
IF (JM1l .LT. 1) GO TO 40 


COMPUTE U(1,J) ,-T=1, ssc,a-1 
me 35 I=-1,JM1 
SUM = LU(I,J) 
IM1 = I-1 


AI = 


mw(EDGT -EO. 0) GO TO 2S 


WITH ACCURACY TEST 
DABS (SUM) 


jig 


LUDAO0630 
LUDAO 64 0 
LUDAO650 
LUDA0660 
LUDAO 670 
LUDAO06 80 
LUDA0690 
LUDAO700 
LUDA0710 
LUDAO7 20 
LUDAO0730 
LUDAO 740 
LUDAO750 
LUDAO760 
LUDAO 770 
LUDAO 780 
LUDA0790 
LUDAO 800 
LUDAO810 
LUDAO8 20 
LUDAO 830 
LUDAO 840 
LUDAO 850 


LUDAO0910 
LUDA0920 
LUDA0930 
LUDA0940 
LUDAO0950 
LUDAO0960 
LUDAO9 70 
LUDAOY 80 
LUDA099 0 
LUDA1000 
LUDA1010 
LUDA1020 
LUDA1030 
LUDA1040 
LUDA1050 
LUDA1060 
LUDA1070 
LUDA10 80 
LUDA1090 
LUDA1100 
LUDA1110 
LUDA1120 
LUDA1130 
LUDA1140 
LUDA1150 
LUDA1160 
LUDA1170 
LUDA1180 
LUDA1190 
LUDA1200 
LUDA1210 
LUDA1220 
LUDA12 30 
LUDA1 240 


i 


20 


25 


30 


22 
40 


45 


5 


SS 


60 


65 


70 


fi: 


WI = ZERO 
IF (IMs Lie.) cent oma0 
DO 15 K=1,IM1 

T = LU(I,K) *LU(K,J) 

SUM = SUM-T 

WI = WI+DABS(T) 
CONTINUE 
LU(I,J) = SUM 
WI = WI+DABS (SUM) 
IF (AI .EQ. ZERO) AI = BIGA 
TEST = WI/AI 
IF (TEST .GT. WREL) WREL = TEST 
GO TO 35 

WITHOUT ACCURACY 

IF (IM1 .LiS 1)8G®8 Toss 
DO 30 K=1,IM1 

SUM = SUM-LU(I,K) *LU(K,J) 
CONTINUE 


LOX) = or 


CONTINUE 


P 


DO 


= ZERO 
COMPUTE U (UigeeeAND L(1,J), Ll=J+1, 28 


70 I=J,N 
SUM = LU(I,J) 
IF (IDGT .EQO. 0) GO TO 5&5 
WITH ACCURACY TEST 

DABS (SUM) 
ZERO 
IF (JM1 .LT. 1) GO TO 50 
DO 45 K=1,JM1 

T = LU(L,K) *LU(kg) 

SUM = SUM-T 

WI = WI+DABS(T) 
CONTINUE 
LU(I,J) = SUM 
WI = WI+DABS (SUM) 
IF (AI .EQ. ZERO) AI = BIGA 
TEST = WI/AI 
IF (TEST .GT. WREL) WREL = TEST 
GO TO 65 


by 
H 
nou 


WITHOUT ACCURACY TEST 

IF (JM1l .LT. TeGOrTe7cs 
DO 60 K=1,JM1 

SUM = SUM-LU(1i7) iu) 
CONTINUE 
LU(i,u) = SUM 
O = EQUIL (I) *DABS (SUM) 
IF (P .GE. Q) GO TO 70 
Pe=rO 
IMAX = I 


CONTINUE 


IF 
IF 


D1 
DO 


TEST FOR ALGORITHMIC SINGULARITY 
(RN+P .EQ. RN) GO TO 110 
(J .EQ. IMAX) GO TO 80 
INTERCHANGE ROWS J AND IMAX 
= = Di 
715 Kel, N 
P = LU(IMAX, K) 
LU (IMAX,K) = LU(d,K) 
LU (J, K) = Pp 


CONTINUE 


‘EQUIL(IMAX) = EQUIL(J) 
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LUDA1250 
LUDA1260 
LUDA1270 
LUDA1280 
LUDA1290 
LUDA1300 
LUDA1310 
LUDA1320 
LUDA1330 
LUDA1340 
LUDA1350 
LUDA1360 
LUDA1370 
LUDA13 80 
LUDA1390 
LUDA1400 
LUDA1410 
LUDA1420 
LUDA14 30 
LUDA1440 
LUDA1450 
LUDA1460 
LUDA1470 
LUDA14 80 
LUDA1490 
LUDA1500 
LUDA1510 
LUDA1520 
LUDA1530 
LUDA1540 
LUDA1550 
LUDA1560 
LUDA1570 
LUDA15 80 
LUDA1590 
LUDA1600 
LUDA1610 
LUDA1620 
LUDA1630 
LUDA1640 
LUDA1650 
LUDA1660 
LUDA1670 
LUDA1680 
LUDA1690 
LUDA1700 
LUDA1710 
LUDA1720 
LUDA1730 
LUDA1740 
LUDA1750 
LUDA1760 
LUDA1770 
LUDA17 80 
LUDA1790 
LUDA1 800 
LUDA1810 
LUDA1 820 
LUDA1830 
LUDA1 840 
LUDA1850 
LUDA1860 


a Oe eee eee eee ee_rlee 


———————E 


80 


a5 


90 


a5 


100 
£OS 


210 


9000 


9005 


IPVT(J) = IMAX 
Dl = D1*LU(J,Jd) 
IF (DABS(D1) .LE. 
Dl = D1*SIXTH 
D2 = D2+FOUR 
GO TO 85 
IF (DABS(D1) .GE. 
Dl = D1*SIXTN 
D2 = D2-FOUR 
GO TO 90 
CONTINUE 
JP1l = J+l 
IF (JP1 .GT. N) GO TO 105 

DIVIDE BY PIVOT ELEMENT U(J,J) 


ONE) GO TO 90 


SIXTH) GOMBGg25 


P = LU(J,J) 
DO 100 I=JP1,N 
LU(I,J) = LU(I,J)/P 

CONTINUE 

CONTINUE 
PERFORM ACCURACY TEST 

IF (IDGT .EQ. 0) GO TO 9005 
P = 3*N+3 
WA = P*WREL 
IF (WA+10.D0**(-IDGT) .NE. WA) GO TO 9005 


IER = 34 
GO TO 9000 
ALGORITHMIC SINGULARITY 
fae = 129 
D1 = ZERO 
D2 = ZERO 
CONTINUE 


PRINT ERROR 
CALL UERTST (IER, 6HLUDATF) 
RETURN 
END 
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LUDA1870 
LUDA18 80 
LUDA1890 
LUDA1I 00 
LUDA1910 
LUDA19 20 
LUDA1930 
LUDA1940 
LUDA1950 
LUDA1960 
LUDA1970 
LUDA19 80 
LUDA1990 
LUDA2000 
LUDA2 010 
LUDA2020 
LUDA2030 
LUDA2040 
LUDA2050 
LUDA2 060 
LUDA2 070 
LUDA2 080 
LUDA2 090 
LUDA2100 
LUDA2110 
LUDA2120 
LUDA2130 
LUDA2140 
LUDA2150 
LUDA2160 
LUDA2170 
LUDA2180 
LUDA2190 
LUDA2200 
LUDA2210 
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IMSL ROUTINE NAME - LUELMF LUEFO010 


LUEFO0020 

wm SR nt et ir Bi en i aa ee LUEF0030 
LUEF0040 

COMPUTER - IBM/DOUBLE LUEFO050 
LUEFO060 

LATEST REVISION - JANUARY 1, 1978 LUEF0070 
LUEFO0080 

PURPOSE -~ ELIMINATION PART OF SOLUTION OF AX=B LUEFO090 
(FULL STORAGE MODE) LUEFO100 

LUEFO110 

USAGE - CALL LUELMF (A,B, IPVT,N, IA, X) LUEF0O120 
LUEFO130 

ARGUMENTS A - A = LU (THE RESULT COMPUTED IN THE IMSL LUEFO0140 
ROUTINE LUDATF) WHERE L IS A LOWER LUEFO150 

TRIANGULAR MATRIX WITH ONES ON THE MAIN LUEF0160 

DIAGONAL. U IS UPPER TRIANGULAR. L AND U LUEF0O170 

ARE STORED AS A SINGLE MATRIX A AND THE LUEFO180 

UNIT DIAGONAL OF L IS NOT STORED. (INPUT) LUEFO190 

B - B IS A VECTOR OF LENGTH N ON THE RIGHT HAND LUEFO0200 

SIDE OF THE EQUATION AX=B. (INPUT) LUEFO210 

IPVT - THE PERMUTATION MATRIX RETURNED FROM THE LUEFO220 

IMSL ROUTINE LUDATF, STORED AS AN N LENGTH LUEFO230 

VECTOR. (INPUT) LUEFO240 

N - ORDER OF A AND NUMBER OF ROWS IN B. (INPUT) LUEF0250 

IA - ROW DIMENSION OF A EXACTLY AS SPECIFIED IN LUEF0260 

THE DIMENSION STATEMENT IN THE CALLING LUEF0270 

PROGRAM. (INPUT) LUEFO0280 

Xx ~ THE RESULT X-. (OURPUT) LUEFO290 

LUEFO0300 

PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 LUEF0310 
- SINGLE /H36,H48,H60 LUEF0320 

LUEFO0330 

REQD. IMSL ROUTINES - NONE REQUIRED LUEFO0340 
LUEFO350 

NOTATION - INFORMATION ON SPECIAL NOTATION AND LUEFO0360 
CONVENTIONS IS AVAILABLE IN THE MANUAL LUEF0370 

INTRODUCTION OR THROUGH IMSL ROUTINE UHELP LUEFO0380 

LUEFO390 

COPYRIGHT - 1978 BY IMSL,” ING? ALL RIGHES RESERVED: LUEF0400 
LUEF0410 

WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN LUEF0420 
APPLIED TO THIS CODE. NO OTHER WARRANTY, LUEF0430 

EXPRESSED OR IMPLIED, IS APPLICABLE. LUEFO0O440 

LUEFO450 

ene aes LUEF0460 
LUEF0470 

SUBROUTINE LUELMF (A,B, IPVT,N,IA,X) LUEF0480 
LUEF0490 

DIMENSION ACTA, 1), Bl) tevin) ex) LUEFO500 
DOUBLE PRECISION A,B,X, SUM LUEFO510 
FIRST EXECUTABLE STATEMENT LUEFO520 

SOLVE LY = B FOR Y LUEF0530 

DOws: f= LUEFO0540 

Sy er) LUEFO550 
IW = 0 LUEFO0560 
DO 20 f-27) LUEF0570 
IP = IP) LUEFO580 

SUM = X(IP) LUEF0590 
X(IP) = X(TI) LUEFO600 

IF (IWe- EO. 0) GOsTGers LUEF0610 

Mi = eed LUEFO0620 


80 


Posto J=I1W,1IM1 
See SUM-A(I, Jd) (J) 
10 CONTINUE 


sO TO 20 
eS Te (SUM .NE. 0.D0) IW = I 
20 X(I) = SUM 


SOLVE UX = Y FOR™Z 
fm 50 IB=1,N 


T = N+1-IB 
TIP1 = I+1l 
SUM = X(I) 


toes. -GT. N) GO TO 30 
co 25 J=2P1,N 
Sil = SUM-A(T,d) *X (J) 
22 CONTINUE 
meertt) = SUM/A(I,I) 
RETURN 
END 


om 


LUEF0630 
LUEFO640 
LUEFO6S0 
LUEFO660 
LUEF0O670 
LUEFO680 
LUEFO690 
LUEFO0700 
LUEFO710 
LUEFO720 
LUEF0730 
LUEFO740 
LUEEOT7S0 
LUEFO760 
LUEF0770 
LUEF0780 
LUEFO790 
LUEFO800 
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IMSL ROUTINE NAME 


COMPUTER 


LATEST REVISION 


PURPOSE 


USAGE 


ARGUMENTS A 


N 


NLC 


IA 


IB 


IJOB 


XL 


IER 


- JANUARY 1, 


- WORK AREA OF DIMENSION N* (NLC+1) . 
NLC*N LOCATIONS OF XL CONTAIN COMPONENTS OF 


LEQT1B 


IBM/DOUBLE 


17 8 


LINEAR EQUATION SOLUTION - BAND STORAGE 


MODE - 


- CALL LEQTIB 


SPACE ECONOMIZER SOLUTION 


(A,N,NLC,NUC,IA,B,M,1IB,IJOB,XL, 


IER) 


INPUT/OUTPUT MATRIX OF DIMENSION N BY 


(NUC+NLC+1). SEE PARAMETER IJOB. 

ORDER OF MATRIX A AND THE NUMBER OF ROWS IN 
B. (INPUT) 

NUMBER OF LOWER CODIAGONALS IN MATRIX A. 
(INPUT) 

NUMBER OF UPPER CODIAGONALS IN MATRIX A. 
(INPUT) 


ROW 


SPECIFIED IN THE DIMENSION STATEMENT IN THE 
CALLING PROGRAM. 


DIMENSION OF MATRIX A EXACTLY AS 


(INPUT) 


INPUT/OUTPUT MATRIX OF DIMENSION N BY M. 


ON INPUT, B CONTAINS THE M RIGHT-HAND SIDES 
OF THE EQUATION AX = 
SOLUTION MATRIX X REPLACES B. 


B 


NUMBER OF RIGHT HAND SIDES 


B. ON OUTPUT, THE 
IF IJOB = 1, 
LS NOt USED. 


(COLUMNS IN B). 


(INPUT) 


ROW 


SPECIFIED IN THE DIMENSION STATEMENT IN THE 
CALLING PROGRAM. 
INPUT OPTION PARAMETER. IJOB = 


Z 


DIMENSION OF MATRIX B EXACTLY AS 
(INPUT) 


= 0, FACTOR THE MATRIX A AND SOLVE THE 
EQUATION AX = 


WHERE A IS ASSUMED TO BE AN N BY N BAND 
MATRIX. A IS STORED IN BAND STORAGE MODE 
AND THEREFORE HAS DIMENSION N BY 
(NLC+NUC+1) . ON OUTPUT, A IS REPLACED 
BY THE U MATRIX OF THE L-U DECOMPOSITION 
OF A ROWWISE PERMUTATION OF MATRIX A. U 
IS STORED IN BAND STORAGE MODE. 

= 1, FACTOR THE MATRIX A. A CONTAINS THE 
SAME INPUT/OUTPUT INFORMATION AS IF 
LI eee 

= 2, SOLVE THE EQUATION AX = B. THIS 
OPTION IMPLIES THAT LEQT1B HAS ALREADY 
BEEN CALLED USING IJOB = O OR 1 SO THAT 
THE MATRIX A HAS ALREADY BEEN FACTORED. 
IN THIS CASE, OUTPUT MATRICES A AND XL 
MUST HAVE BEEN SAVED FOR REUSE IN THE 
CALL TO LEQT1B. 

THESE ERS] 


THE L MATRIX OF THE L-U DECOMPOSITION OF A 
ROWWISE PERMUTATION OF A. THE LAST N 
LOCATIONS CONTAIN THE PIVOT INDICES. 


ERROR PARAMETER. 


(OUTPUT) 


Se 


I IMPLIES WHEN 


B. ON INPUT, A CONTAINS THE 
COEFFICIENT MATRIX OF THE EQUATION AX = B, 


LE1B0010 
LE1B0020 


LE1B0040 
LE1B0050 
LE1B0060 
LE1B0070 
LE1BO0080 
LE1B0090 
LE1B0100 
LESS Oo 
LE1B0120 
LE1B0130 
LE1B0140 
LE1B0150 
LE1B0160 
LE1B0170 
LE1B0180 
LE1B0190 
LE1B0200 
LE1B0210 
LE1B0220 
LE1B0230 
LE1B0240 
LE1B0250 
LE1B0260 
LE1B0270 
LE1B0280 
LE1B0290 
LE1B0300 
LE1B0310 
LE1B0320 
LE1B0330 
LE1B0340 
LE1B0350 
LE1B0360 
LE1B0370 
LE1B0380 
LE1B0390 
LE1B0400 
LE1B0410 
LE1B0420 
LE1B0430 
LE1B0440 
LE1B0450 
LE1B0460 
LE1B0470 
LE1B0480 
LE1B0490 
LE1B0500 
LE1B0510 
LE1B0520 
LE1B0530 
LE1B0540 
LE1B0550 
LE1B0560 
LE1B0570 
LE1B0580 
LE1B0590 
LE1B0600 
LE1B0610 
LE1B0620 
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LS 
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PRECISION/HARDWARE 


TERMINAL ERROR 
IER = 129 INDICATES THAT MATRIX A IS 
ALGORITHMICALLY SINGULAR. (SEE THE 
CHAPTER L PRELUDE). 


- SINGLE AND DOUBLE/H32 
- SINGLE/H36,H48,H60 


fae. IMSL ROUTINES - UERTST,UGETIO 

NOTATION - INFORMATION ON SPECIAL NOTATION AND 
CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 
COPYRIGHT ~el7g?o Bye lM>ob, onc. Sob RIGHTS RESERVED. 


WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 


EXPRESSED OR IMPLIED, IS APPLICABLE. 


SeeROULINE LEOQTIB (A,N,NLC,NUC,1IA,B,M, IB, IJOB, XL, IER) 


DIMENSION Peek 1) xL (N, 1) 28 (0B, 1) 
DOUBLE PRECISION A,XL,B,P,Q,ZERO,ONE, RN 
BATA ZERO/0.D0/,ONE/1.0D0/ 
FIRST EXECUTABLE STATEMENT 
IER = 0 
JBEG NLC+1 


NLC1 = JBEG 
me (1JOB .EQ. 2) GO TO 80 


pay = N 
RESTRUCTURE THE MATRIX 
FIND RECIPROCAL OF THE LARGEST 
ABSOLUTE VALUE IN ROW I 

oo | 

NC = JBEG+NUC 

NN = NC 

JEND = NC 

fen .E0O. 1 .OR. NLC .EQ. 0) GO TO 25 

fs 1 

P = ZERO 

DO 10 J = JBEG, JEND : 


me, K) = A(I,J) 
Q = DABS(A(I,K)) 
IF (0 .GT. P) P=Q 
K = K+l 
CONTINUE 
IF (P .EQ. ZERO) GO TO 135 
XL(I,NLC1) = ONE/P 
IF (K .GT. NC) GO TO 20 
DO 15 J = K,NC 
A(I,J) = ZERO 
CONTINUE 
I = I+l1 
JBEG = JBEG-1 
IF (JEND-JBEG .EQ. N) JEND = JEND-1 
IF (I .LE. NLC) GO TO 5 
JBEG = I 
NN = JEND 
JEND = N-NUC 


oo 


LE1B0630 
LE1B0640 
LE1B0650 
LE1B0660 
LE1B0670 
LE1B0680 
LE1B0690 
LE1B0700 
LE1B0710 
LE1B0720 
LE1B0730 
LE1B0740 
LE1B0750 
LE1B0760 
LE1B0770 
LE1B0780 
ETE G72 0 
LE1B0800 
LE1B0810 
LE1B0820 


we en ee ee ee ee ee ee ee LE1B0830 


LE1B0840 
LE1B0850 
LE1B0860 
LE1B0870 
LE1B08 80 
LE1B0890 
LE1B0900 
LE1B0910 
LE1B0920 
LE1B0930 
LE1B0940 
LE1B0950 
LE1B0960 
LE1B0970 
LE1B0980 
LE1B0990 
LE1B1000 
LE1B1010 
LE1B1020 
LE1B1030 
LE1B1040 
LE1B1050 
LE1B1060 
LE1B1070 
LE1B1080 
LE1B1090 
LE1B1100 
LE1B1110 
LE1B1120 
LEtBt130 
LE1B1140 
LE1B1150 
LE1B1160 
LE1B1170 
LE1B1180 
LE1B1190 
LEIBLZ200 
LE1B1210 
LE1B1220 
LE1B1230 
LE1B1240 


DO 40 I = JBEG,N LEDEYZS¢ 





P = ZERO LE1B1260 
DO 30 J = 1,NN LE1B1270 
Q = DABS(A(I,J)) LE1B1280 
IF (Q .GT. P) P= Q LE1B1290 
30 CONTINUE LE1B1300 
IF (P .EQ. ZERO) GO TO 135 LE1B1310 
XL(I,NLC1) = ONE/P LE1B1320 
IF (I .EQ. JEND) GO TO 37 LE1B1330 
IF (I .LT. JEND) GoeTro 40 LE1B1340 
K = NN+l LE1B1350 
DO 35 J = K,NC LE1B1360 
A(I,d) = ZERO LE1B1370 
35 CONTINUE LE1B1380 
37 NN = NN-1 LE1B1390 
40 CONTINUE LE1B1400 
L = NLC LE1B1410 
L-U DECOMPOSITION LE1B1420 
DO 75 K = 1,N LE1B1430 
P = DABS(A(K,1))*XL(K,NLC1) LE1B1440 
Ir=K LE1B1450 
IF (L .LT. N) L = L+l1 LE1B1460 
Kl = K+l LE1B1470 
IF (K1 .GT. L) GO TO 50 LE1B1480 
DO 45 J = K1,L LE1B1490 
Q = DABS (A(J,1))*XL(J,NLC1) LE1B1500 
IF (Q .LE. P) GO TO 45 LE1B1510 | 
P=0Q LE1B1520 . 
I=eJ LE1B1530 ) 
45 CONTINUE LE1B1540 
SO XL(I,NLC1) = XIV¢he timed) LE1B1550 
XL(K,NUGi) = 1 LE1B1560 
SINGULARITY FOUND LE1B1570 
Q = RN+P LE1B1580 
IF (Q .EQ. RN) GO TO 135 LE1B1590 
INTERCHANGE ROWS I AND K LE1B1600 
IF (K .EQ. I) GO TO 60 LE1B1610 
DO 55 J = 1,NC LE1B1620 
P = A(K,J) LE1B1630 
RK gd) = oe LE1B1640 
A(I,J) = P LE1B1650 
55 CONTINUE LE1B1660 
60 IF (Kl .GT. L} GO Toms LE1B1670 
BO 7Or io =chler, : LE1B1680 
P= A(I,1)/A(kX, LE1B1690 
IK = I-K LE1B1700 
XL(K1,IK) = P LE1B1710 | 
DO 65 J = 2,NC LE1B1720 
A(I,J-1) = A( Pyar alke os) LE1B1730 
65 CONTINUE LE1B1740 
A(I,NC) = ZERO LE1B1750 | 
70 CONTINUE LE1B1760 
75 CONTINUE LE1B1770 
IF (IJOB .EQ. 1) GO TO 9005 LE1B1780 
FORWARD SUBSTITUTION LE1B1790 
80 L = NLC LE1B1800 
DO 105 K = 1,N LE1B1810 
T = XEtK NLC) LE1B1820 
IF (I .EO. K) GO Ten90 LE1B1830 
DO 85 J =1,M LE1B1840 
P = B(K,d) LE1B1850 
B(K,J) = B(I,Jd) LE1B1860 


84 


SE 


| 
i 
i 
i 
Hi 
4 
; 





90 


95 
100 
EOD 


10 
aS 


12,0 
125 


22 
9000 


9005 


Bel =) P 
CONTINUE 
mee ee, N) Gb = L+l 
Kl = K+l 
meen? .GT. L) GO TO 105 
emo I = K1 ,L 
eo =m - Kk 
P = XL(K1, 1K) 
ae Ss J =" bam 
Bal, J) = B(1,J) -P*B(K,0) 
CONTINUE 
CONTINUE 
CONTINUE 
BACKWARD SUBSTITUTION 
JBEG = NUC+NLC 
mm 125 J = 1,M 


L=l1 

Kl = N+tl 

DO 120 I =1,N 
K = K1-I 
P = B(K,J) 


mw (L .EO.*1) “GO TO 115 
ie Teo KK = 2,L 
IK = KK+K 
P = P-A(K,KK) *B(IK-1,J) 
CONTINUE 
Bikeway (kK, 2) 
IF (L .LE. JBEG) L = L+tl 
CONTINUE 
CONTINUE 
GO TO 9005 
IER = 129 
CONTINUE 
CALL UERTST (IER, 6HLEQT1B) 
RETURN 
END 


oo 


LE1B1870 
LE1B1880 
LE1B1890 
LE1B1900 
Eee 2 10 
LE1B1920 
Eee .930 
LE1B1940 
LE1B1950 
LE1B1960 
LE1B1970 
LE1B1980 
LE1B1990 
LE1B2000 
LE1B2010 
LE1B2020 
LE1B2030 
LE1B2040 
LE1B2050 
LE1B2060 
LE1B2070 
LE1B2080 
LE1B2090 
LE1B2100 
LEYE2110 
LE1B2120 
LE1B2130 
LE1B2140 
LE1B2150 
LE1B2160 
LE1B2170 
LE1B2180 
DELB2Z1I0 
LE1B2200 
LE1B2210 
bELTBZ220 
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IMSL ROUTINE NAME - UERTST 
COMPUTER - IBM/SINGLE 
LATEST REVISION - JUNE 1, 1982 
PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION 
USAGE - CALL UERTST (IER, NAME) 
ARGUMENTS IER - ERROR PARAMETER. (INPUT) 
IER = I+J WHERE 
I = 128 IMPLIES TERMINAL ERROR MESSAGE, 
I = 64 IMPLIES WARNING WITH FIX MESSAGE, 
I = 32 IMPLIES WARNING MESSAGE. 
J = ERROR CODE RELEVANT TO CALLING 
ROUTINE. 
NAME - A CHARACTER STRING OF LENGTH SIX PROVIDING 
THE NAME OF THE CALLING ROUTINE. (INPUT) 
PRECISION/HARDWARE - SINGLE/ALL 
REQD. IMSL ROUTINES - UGETIO, USPKD 
NOTATION - INFORMATION ON SPECIAL NOTATION AND 
CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 
REMARKS THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN 
TO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT 
NUMBER CAN BE DETERMINED BY CALLING UGETIO AS 
FOLLOWS. . CALL UGETIO(1,NIN,NOUT) . 
THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING 
UGETIO AS FOLLOWS... 
NIN = 0 
NOUT = NEW OUTPUT UNIT NUMBER 
CALL UGETIO(3,NIN,NOUT) 
SEE THE UGETIO DOCUMENT FOR MORE DETAILS. 
COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. 
WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 
SUBROUTINE UERTST (IER, NAME) 
SPECIFICATIONS FOR ARGUMENTS 
INTEGER IER 
INTEGER NAME (1) 
SPECIFICATIONS FOR LOCAL VARIABLES 
INTEGER I, IEQ, IEQDF, IOUNIT, LEVEL, LEVOLD, NAMEQ (6), 
* NAMSET (6) , NAMUPK (6) , NIN, NMTB 
DATA NAMSET/1HU,1HE,1HR,1HS,1HE,1HT/ 
DATA NAMEQ/6*1H / 
DATA LEVEL/4/,IEQDF/0/,IEQ/1H=/ 
UNPACK NAME INTO NAMUPK 
FIRST EXECUTABLE STATEMENT 
CALL USPKD (NAME, 6,NAMUPK,NMTB) 


86 


UERTO010 
UVERTO020 


UERTOOSO 
UERTOO60 
UERTO070 
UERTOO80 
UERTOO090 
UERTO100 
UERTO110 
UERTO120 
UERTO130 
UERTO140 
UERTO150 
UERTO160 
UERTO170 
UERTO180 
UERTO190 
UERTO200 
UERTO210 
UERTO220 
UERTO230 
UERTO240 
UERTO250 
UERT0260 
UERT0270 
UERTO280 
UERTO290 
UERT0300 
UERTO310 
UERTO320 
UERTO330 
UERT0340 
UERTO350 
UERT0360 
UERT0370 
UERTO3 80 
UERT0390 
UERT0O400 
UERT0410 
UERT0420 
UERTO430 
UERT0440 
UERTO450 
UERT0460 
UERT0470 


UERTOS00 
UERTO510 
UERTO520 
UERTO530 
UERTO540 
UERTOS5S50 
UERTOS560 
UERTO570 
UERTO580 
UERTOS90 
UERTO600 
UERTO610 
UERTO620 
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H 
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Cerra 


co 


iS 


20 


a 


20 
35 
40 
45 


210 


ao 


60 
65 


GET OUTPUT UNIT NUMBER 
CALL UGETIO(1,NIN, IOUNIT) 
CHECK IER 
IF (IER.GT.999) GO TO 25 
IF (IER.LT.-32) GO TO 55 
IF (IER.LE.128) GO TO 5 
me (LEVEL.LT.1) GO TO 30 
PRINT TERMINAL MESSAGE 
IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER, NAMEQ, IEQ, NAMUPK 
IF (IEQDF.EQ.0) WRITE (IOUNIT,35) IER, NAMUPK 
GO TO 30 
IF (IER.LE.64) GO TO 10 
ma (LEVEL.LT.2) GO TO 30 
PRINT WARNING WITH FIX MESSAGE 
IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ, IEQ, NAMUPK 
IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAMUPK 
GO TO 30 
IF (IER.LE.32) GO TO 15 
PRINT WARNING MESSAGE 
me (LEVEL.LT.3) GO TO 30 
IF (IEQDF.EQ.1) WRITE (IOUNIT,45) IER,NAMEQ, IEQ, NAMUPK 
IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAMUPK 
GO TO 30 
CONTINUE 
CHECK FOR UERSET CALL 
DO 20 I=1,6 
IF (NAMUPK(I) .NE.NAMSET(I)) GO TO 25 
CONTINUE 
LEVOLD = LEVEL 
LEVEL = IER 
IER = LEVOLD 


IF (LEVEL.LT.0) LEVEL = 4 
IF (LEVEL.GT.4) LEVEL = 4 
GO TO 30 
CONTINUE 


me (LEVEL.LT.4) GO TO 30 


PRINT NON-DEFINED MESSAGE 
IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER, NAMEQ, IEQ, NAMUPK 


IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER, NAMUPK 
mODF = 0 
RETURN 


PORMAT (19H *** TERMINAL ERROR,10X,7H(IER = ,13, 
20H) FROM IMSL ROUTINE ,6A1,Al1,6A1) 
FORMAT (27H *** WARNING WITH FIX ERROR,2X,7H(IER = ,13, 
20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 
FORMAT (18H *** WARNING ERROR,11X,7H(IER = ,13, 
20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 
FORMAT (20H *** UNDEFINED ERROR,9X,7H(IER = ,I15, 
20H) FROM IMSL ROUTINE ,6A1,Al1,6A1) 


SAVE P FOR P = R CASE 

P IS THE PAGE NAMUPK 

R IS THE ROUTINE NAMUPK 
IEQDF = 1 
DO 60 I=1,6 
NAMEQ(I) = NAMUPK(I) 
RETURN 
END 


oy 


UERTO630 
UERT0640 
UERTO6S50 
JERTO660 
UERT0670 
UERTO680 
UERTO690 
UERTO700 
UERTO710 
UERT0O720 
UERTO730 
UERTO740 
UERTO750 
UERTO760 
UERT0770 
UERTO780 
UERTO790 
UERTO800 
UERTO810 
UERTO820 
UERTO830 
UERTO840 
UERTO850 
UERTO860 
UERTO870 
UERTO880 
UERTO890 
UERTO900 
UERTO910 
UERTO920 
UERTO930 
UERTO940 
UERTO950 
UERTO960 
UERT0O970 
UERTO980 
UERTO990 
UERT1000 
UERT1010 
UERT1020 
UERT1030 
UERT1040 
UERT1050 
UERT1060 
UERT1070 
UERT1080 
UERT1090 
UERT1100 
UERT1110 
UERT1120 
WERT lus 0 
UERT1140 
UERT1150 
UERT1160 
UERTI1T7O 
UERT1180 
UVERT1190 
UERT1200 
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IMSL ROUTINE NAME 


COMPUTER 


LATEST REVISION 


PURPOSE 
USAGE 
ARGUMENTS IOPT 
NIN 
NOUT 
PRECISION/HARDWARE 
REQD. IMSL ROUTINES 
NOTATION 
REMARKS 
IOPT=3, 
COPYRIGHT 
WARRANTY 


- UGETIO 
- IBM/SINGLE 
- JUNE 1, 29e1 


- TO RETRIEVE CURRENT VALUES AND TO SET NEW 
VALUES FOR INPUT AND OUTPUT UNIT 
IDENTIFIERS. 


- CALL UGETIO(IOPT, NIN, NOUT) 


- OPTION PARAMETER. 
IE GOer i 


(INPUT) 
THE CURRENT INPUT AND OUTPUT 


UNIT IDENTIFIER VALUES ARE RETURNED IN NIN 
AND NOUT, RESPECTIVELY. 

IF IOPT=2, THE INTERNAL VALUE OF NIN IS 
RESET FOR SUBSEQUENT USE. 


Lee LOPL=3, 


THE INTERNAL VALUE OF NOUT IS 


RESET FOR SUBSEQUENT USE. 


INPUT UNIT IDENTIFIER. 
OUEPULT Ir TORR —i. 


INPUT IF IOPT=2. 


- OUTPUT UNIT IDENTIFIER. 


OUTEUT fe ort i, 


TNPUL IF IOP 22 


- SINGLE /ALL 


- NONE REQUIRED 


INFORMATION ON SPECIAL NOTATION AND 


CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 


EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT 
OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT 


IDENTIFIER VALUES. 


IF UGETIO IS CALLED WITH IOPT=2 OR 


NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED. 
SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. 


1978 BY IMSL, INC. ALL RIGHTS RESERVED. 


IMSL WARRANTS ONLY THAT 
APPLIED TO THIS CODE. 
EXPRESSED OR IMPLIED, 


IMSL TESTING HAS BEEN 
NO OTHER WARRANTY, 
IS APPLICABLE. 


SUBROUTINE UGETIO(IOPT, NIN, NOUT) 


INTEGE 


R 


INTEGER 


DATA 


ask 
IF 
rE 
NEN — 
NOUT = 
GO TO 


NIND 
NOUTD 
9005 


SPECIFICATIONS FOR ARGUMENTS 


heer, NIN, NOUT 


SPECIFICATIONS FOR LOCAL VARIABLES 


NIND, NOUTD 
NIND/5/,NOUTD/6/ 


(IOPT.EQ.3) GO TO 10 
(IOPT.EQ.2) GO TO 5 
(IOPT .NE.1) GO Tes 9005 


FIRST EXECUTABLE STATEMENT 


88 


UGETO010 
UGET0020 


wwe ee nn sn oo wi ie UGET0030 


UGETO040 
UGETOOS50 
UGETOO60 
UGET0070 
UGETO080 
UGETO090 
UGETO100 
UGELOLES 
UGETO120 
UGETOTR® 
UGETO140 
UGETO150 
UGETO160 
UGET0170 
UGETO180 
UGETO190 
UGETO200 
UGETO210 
UGET0O220 
UGETO230 
UGET0240 
UGETO250 
UGET0260 
UGET0270 
UGET0280 
UGET0290 
UGET0300 
UGET0310 
UGET0320 
UGET0330 
UGET0340 
UGETO350 
UGET0360 
UGET0370 
UGET03 80 
UGET0390 
UGET0400 
UGET0410 
UGET0420 
UGET0430 
UGET0440 
UGETO450 
UGET0460 
UGET0470 


ee UGET04 80 


UGET0490 
UGETO500 
UGETOS510 
UGETO520 
UGET0530 
UGET0540 
UGETOS550 
UGETO560 
UGET0570 
UGETO580 
UGETO590 
UGETO0600 
UGETO610 
UGET0620 





89 


UGET0630 
UGET0640 
UGET0650 
UGET0660 
UGET0670 
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USPK0010 
USPK0020 


USPK0040 
USPKO0050 
USPKO0060 
USPK0070 
USPK0080 
USPKO090 
USPK0100 
USPK0110 
USPK0120 
USPK0130 
USPK0140 
USPK0150 
USPK0160 
USPK0170 
USPK0180 
USPK0190 
USPK0200 
USPK0210 


IMSL ROUTINE NAME - USPKD 
wwe ee ee ee ee eee eee ee e+ + ------------ USPK0030 
COMPUTER - IBM/SINGLE 
LATEST REVISION - NOVEMBER 1, 1984 
PURPOSE - NUCLEUS CALLED BY IMSL ROUTINES THAT HAVE 
CHARACTER STRING ARGUMENTS 
USAGE - CALL USPKD (PACKED, NCHARS, UNPAKD , NCHMTB) 
ARGUMENTS | PACKED - CHARACTER STRING TO BE UNPACKED. (INPUT) 
NCHARS - LENGTH OF PACKED. (INPUT) SEE REMARKS. 
UNPAKD - INTEGER ARRAY TO RECEIVE THE UNPACKED 
REPRESENTATION OF THE STRING. (OUTPUT) 
. NCHMTB - NCHARS MINUS TRAILING BLANKS. (OUTPUT) 
PRECISION/HARDWARE - SINGLE/ALL 
REQD. IMSL ROUTINES - NONE 


REMARKS 1. 
IN (Al) 


2. UP TO 129 CHARACTERS MAY BE USED. 


USPKD UNPACKS A CHARACTER STRING INTO AN INTEGER ARRAY 


FORMAT . 
ANY IN EXCESS OF 


THAT ARE IGNORED. 


COPYRIGHT. 


WARRANTY 


- 1984 BY IMSL, INC. ALL RIGHTS RESERVED. 


APPLIED TO THIS CODE. 
EXPRESSED OR IMPLIED, 


NO OTHER WARRANTY, 
IS APPLICABLE. 


IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 


USPKO0220 
USPK0230 
USPK0240 
USPK0250 
USPK0260 
USPK0270 
USPK0280 
USPK02390 
USPK0300 
USPK0310 
USPK0320 
USPK0330 
USPK0340 


mw ww ein mmm a a a i USPK0350 


SUBROUTINE USPKD 


(PACKED, NCHARS , UNPAKD , NCHMTB) 


SPECIFICATIONS FOR ARGUMENTS 
INTEGER NC, NCHARS , NCHMTB 
LOGICAL*1 UNPAKD (1) , PACKED (1) , LBYTE, LBLANK 
INTEGER*2 IBYTE, IBLANK 
EQUIVALENCE (LBYTE, IBYTE) 
DATA LBLANK /1H / 
DATA IBYTE /1H / 
DATA IBLANK /1H / 
INITIALIZE NCHMTB 
NCHMTB = 0 
RETURN IF NCHARS IS LE ZERO 
IF (NCHARS .LE.0) RETURN 
SET NC=NUMBER OF CHARS TO BE DECODED 
NC = MINO (129,NCHARS) 
NWORDS = NC*4 
J = 1 
DO 110 I = 1,NWORDS,4 
UNPAKD (I) = PACKED (J) 
UNPAKD (I+1) = LBLANK 
UNPAKD (I+2) = LBLANK 
UNPAKD (I+3) = LBLANK 
110 J = J#+l 


CHECK UNPAKD ARRAY AND SET NCHMTB 
BASED ON TRAILING BLANKS FOUND 


DO 200 N = 1,NWORDS,4 
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USPK0360 
USPK0370 
USPK0380 
USPK0390 
USPK0400 
USPK0410 
USPK0420 
USPK0430 
USPK0440 
USPK0450 
USPK0460 
USPK0470 
USPK0480 
USPK0490 
USPKO0500 
USPK0510 
USPK0520 
USPK0530 
USPK0540 
USPK0550 
USPK0560 
USPK0570 
USPK0580 
USPK0590 
USPK0600 
USPK0610 
USPK0620 


as 


nee a tle 





NN = NWORDS - N - 2 
LBYTE = UNPAKD (NN) 
IF(IBYTE .NE. IBLANK) GO TO 210 
200 CONTINUE 
NN = O 
210 NCHMTB = (NN + 3) / 4 
RETURN 
END 


OL 


USPK0630 
USPK0640 
USPKO650 
USPKO660 
USPK0670 
USPK0680 
USPK0690 
USPK0700 


LIST OF REFERENCES 
1. von Braun, W. and Ordway, F., History of Rocketry and Space Travel,3d ed., T. Crowell 
Co., 1974. 
2. Sutton, G.P., Rocket Propulsion Elements, 4th ed., John Wiley & Sons, 1976. 
3. Jahn, R.G., Physics of Electric Propulsion, McGraw-Hill, 1968. 
4. Sutton, G.P., Rocket Propulsion Elements, 5th ed., John Wiley & Sons, 1986. 
5. Stuhlinger, E., Jon Propulsion for Space Flight, McGraw-Hill, 1964. 


6. Myers, R.M., Mantenieks, M.A., and LaPointe, M.R., MPD Thruster Technology, 
NASA-TM-105242, AIAA-91-3568, Sept., 1968. 


7. Myers, R.M., Applied-Field MPD Thruster Geometry Effects, AUAA-91-2342, June, 1991. 


8. Saber, A.J., “Anode Power in the Quasi-Steady MPD (Magnetoplasmadynamic) Thruster," 
Ph.D. Thesis, Princeton Univ., NJ, 1974. 


9. Shih, K.T., Pfender, E., Ibele, W.E., and Eckert, E.R.G., "Experimental Anode Heat- 
Transfer Studies in a Coaxial Arc Configuration," AJAA Journal, V. 6, No. 8, pp.1482- 
1487, August, 1968. 


10. Sanders, N.A., and Pfender, E., "Measurement of the Anode Falls and Anode Heat 
Transfer in Atmospheric Pressure, High Intensity Arcs," J. Appl. Phys., V. 55, No. 3, pp. 
714-722, Feb. ,1984. 


11. Vainberg, L.I., Lyubimov, G.A., and Smolin, G.G., "High Current Discharge Effects 
and Anode Damage in an End-Fire Plasma Accelerator," Sov. Physics, Tech. Phys., V. 23, 
No. 4, pp. 439-443, April,1978. 


12. Hugel, H., "Effects of Self-Magnetic Forces on the Anode Mechanism of a High Current 
Discharge," JEEE Trans. Plasma Sci., V. PS-8, No. 4, pp.437-442, Dec. 1980. 


13. Gallimore, A.D., Kelly, A.J., and Jahn, R.G., "Anode Power Deposition in Quasisteady 
MPD Thrusters,"J. Propulsion & Pwr., V. 8, No. 6, pp. 1224-1231, Dec., 1992. 


pe 


14. Dolson, R.C., and Biblarz, O., "Analysis of the Voltage Drop Arising from a Collision- 
dominated Sheath," J. Appl. Phys., V. 47, No. 12, pp. 5280-5287, Dec., 1976. 


15. Myers, R.M., "Energy Deposition in Low Power Coaxial Plasma Thrusters," Ph.D. 
Thesis, Princeton Univ., NJ, June, 1989. 


16. Oberth, R.C., and Jahn, R.G., “Anode Phenomena in High-Current Accelerators," AIAA 
Journal, V. 10, No. 1, pp. 86-91, Jan., 1972. 


17. Sleziona, P.C., Auweter-Kurtz, M., and Schrade, H.O., "Numerical Codes for 
Cylindrical MPD Thrusters," IEPC 88-038, 20th Int’l. Elec. Prop. Conf., Garmisch- 
Partenkirchen, W. Germany, Oct., 1988. 


18. Sleziona, P.C., Auweter-Kurtz, M., and Schrade, H.O., "Numerical Evaluation of MPD 
Thrusters," AIAA 90-2602, July, 1990. 


19. LaPointe, M.R., "Numerical Simulation of Self-Field MPD Thrusters," AIAA-91-2341, 
NASA-CR-187168, July, 1991. 


20. Subramaniam, V.V., and Lawless, J.L., "Thermal Instabilities of the Anode in a 
Magnetoplamsadynamic Thruster," J. Propulsion & Pwr., V. 6, No. 2, pp. 221-224, Mar., 
1990. 


21. Biblarz, O., "Approximate Sheath Solutions for a Planar Plasma Anode,” JEEE Trans. 
Plasma Sci., V. 19, No. 6, pp. 1235-1243, Dec., 1991. 


22. Biblarz, O., Dolson, R.G., and Shorb, R.C., "Anode Phenomena a Collision-dominated 
Plasma,” J. Appl. Phys., V. 46, No. 8, pp. 3342-3346, Aug., 1975. 


23. Rudolph,L.K. and Pawlik, E.V.,"The MPD Thruster Development Program," AIAA 
Technical Paper 79-2050,in Progress in Aeronautics and Astronautics, Vol. 79, Amer. Inst 
of Aer. & Astro., 1981. 


24. Cann, G.L., and Marlotte, G.L.,"Hall Current Plasma Accelerator," AZAA Journal, V. 
2, No. 7, pp. 1234-1241, Jul., 1964. 


25. Brophy, J., Stationary Plasma Thruster Evaluation in Russia, Summary Report, Jet 
Propulsion Laboratory (JPL) Publication 92-4, March 15, 1992. 


26. Mitchner, M., and Kruger, C.H., Partially Ionized Gases, pp. 128-134, Wiley, 1973. 


27. Cobine, J.D.,Gaseous Conductors, Theory and Engineering Applications, McGraw-Hill, 
1941. 


93 


28. von Engel, A., Jonized Gases, Clarendon Press, 1965. 


29. Chen, F.F., Introduction to Plasma Physics and Controlled Fusion, 2nd ed., p.10, 
Plenum, 1984. 


30. Lin, S., Resler, E., and Kantrowitz, A., "Electrical Conductivity of Highly Ionized 
Argon Produced by Shock Waves," J. Appl. Phys., V. 26, p. 95, Jan., 1955. 


31. Campbell, A., Plasma Physics and Magnetofluidmechanics, p. 161, McGraw-Hill, 1963. 
32. Nasser, E., Fundamentals of Gaseous Ionization and Plasma Electronics, p.412, Wiley, 
1oF le 


33. Ecker, G., "Anode Spot Instability. I. The Homogeneous Short Gap Instability," 
IEEE Trans. Plasma Sct., V. PS-2, No. 3, Sept., 1974. 


34. Biblarz, O. and Riggs, J.F., "Anode Sheath Contributions in Plasma Thrusters," AJAA 
93-2495, Jun., 1993. 


35. Miller, H.C., "Vacuum Arc Anode Phenomena," JEEE Trans. Plasma Sci., V. PS-11, 
No. 2, June, 1983. 


36. Miller, H.C., "Discharge Modes at the Anode of a Vacuum Arc," JEEE Trans. Plasma 
Sci., V. PS-11, No. 3, pp. 122-127, Sep., 1983. 


37. Rich, J.A., Prescott, L.E., and Cobine, J.D., "Anode Phenomena in Metal-Vapor Arcs 
at High Currents," J. Appl. Phys., V. 12, No. 12, pp.587-601, Feb., 1971. 


38. Schuocker, D., "Improved Model for Anode Spot Formation in Vacuum Arcs," JEEE 
Trans. Plasma Sci., V. PS-7, No. 4, pp. 209-216, Dec., 1979. 


39. Gear, C.W., Numerical Initial Value Problems in Ordinary Differential Equations, 
Prentice-Hall, Englewood Cliffs, NJ, 1971. 


40. Burnet, H., Vincent., P., and Rocca Serra, J., "Ionization Mechanism in a Nitrogen 
Glow Discharge," J. Appl. Phys., V. 54, No. 9, pp. 4951-4957, 1983. 


41. Phelps, A.P. and Pitchford, L.C., "Anisotropic Scattering Electrons by N, and its Effect 
on Electron Transport," Phys. Rev. A, V. 31, pp. 2932-2949, 1985. 


42. Myers, R.M., Kelly, A.J.and Jahn, R.G., "Energy Deposition in Low-Power Coaxial 
Thrusters," J. Propulsion, V. 7, No. 5, pp. 732-739, Sep./Oct., 1991. 


94 


ES — EF ——_- —_ — 


43. Biblarz, O., "Thermionic Arc Initiation", 1992 IEEE International Conference on Plasma 
Science, Tampa, FL, June 1992. 


44. Gallimore, A.D., "Anode Power Deposition in Coaxial MPD Thrusters," Ph.D. Thesis, 
Princeton Univ., NJ, Oct., 1992. 


45. Biblarz, O. and Barto, J.L., "Fluid-Dynamic Effects, Including Turbulence, on a High 
Pressure Discharge". Gas Flow and Chemical Lasers, 6th Int. Symposium (S. Rosenwaks, 
Ed.), pp. 34-39, Springer-Verlag, Berlin, 1987. 


46. Park, W. & Choi, D., "Numerical Analysis of MPD Arcs for Plasma Acceleration, " 
IEEE Trans. Plasma Sci., V. PS-15, No. 5, pp. 618-624, Oct., 1987. 


47. Schoeck, P. A., Eckert, E.R.G., and Wutzke, S.A., "An Investigation of the 
Anode Losses in Argon Arcs and their Reduction by Transpiration Cooling," ARL 62-341, 
DTIC AD-278570, April, 1962. 


48. Brady, J.E., General Chemistry, Principles & Structure, 5th ed., Wiley, 1990, p. 223. 


49. Janes, G.S., "Magnetohydrodynamic Propulsion," in Advanced Propulsion Techniques, 
AGARD Proceedings, Aug., 1960, pp. 151-154, Pergamon, 1961. 


50. Connolly, D.J., Sovie, R.J., Michels, C.J., and Burkhart, J.A., "Low Environmental 
Pressure MPD Arc Tests," AIJAA Journal, V. 6, No. 7, pp. 1271-1276, July, 1968. 


51. Subramaniam, V.V., "Fundamental Studies on Erosion in MPD Thrusters," AFOSR 87- 
0360, April, 1992. 


52. Hurwics, H. and Rogan, J. E., "High Temperature Thermal Protection Systems", 
Section 19, , Handbook of Heat Transfer (Rohsenow, W.M., and Hartnett, J.P., Eds) 
McGraw-Hill, 1973. 


53. Kuriki, K. and Suzuki, H., "Quasisteady MPD Arcjet with Anode Gas Injection," AIAA 
79-2058, 14th International Electric Propulsion Conference, Princeton, NJ, Oct., 1979. 


54. Sutton, G.W. and Sherman A., Engineering Magnetohydrodynamics, p. 148, McGraw- 
Hill, 1965. 


55. Choueiri, E.Y., Kelly, A.J., and Jahn, R.G., "Mass Savings Domain of Plasma 
Propulsion for LEO to GEO Transfer," J. Spacecraft and Rockets, V. 30, No. 6, pp. 749- 
754, Nov./Dec., 1993. 


95 


— 


INITIAL DISTRIBUTION LIST 


. Defense Technical Information Center 


Cameron Station 
Alexandria, Virginia 22304-6145 


. Library, Code 52 


Naval Postgraduate School 
Monterey, California 93943-5002 


Chairman 

Department of Aeronautics & Astronautics, Code AA 
Naval Postgraduate School 

Monterey, California 93943-5000 


Professor Oscar Biblarz 

Department of Aeronautics & Astronautics, Code AA/Bi 
Naval Postgraduate School 

Monterey, California 93943-5000 


Professor Fred Schwirzke 
Department of Physics, Code PH/Sw 
Naval Postgraduate School 
Monterey, California 93943-5000 


Commander, Space & Naval Warfare Systems Commamd 
Space Technology Directorate (SPAWAR-40) 

2451 Crystal Drive 

Arlington, Virginia 22245-5200 


. Dr. Roger Myers 


Lewis Research Center 
M.S. SPTD-1 
Cleveland, Ohio 44135-3191 


Dr. James S. Sovey 

SPTD 

NASA Lewis Research Center 
21000 Brookpark Rd. 
Cleveland, Ohio 44135 


96 


10. 


EL . 


2. 


13. 


14. 


ES. 


Dr. Tom Pivirotto 

Jet Propulsion Laboratory 
4800 Oak Grove Dr. 
Pasadena, California 91109 
M.S.125-224 


Dr. John Brophy 

Jet Propulsion Laboratory 
4800 Oak Grove Dr. 
Pasadena, California 91109 
M.S. 125-224 


Dr. Jay Polk 

Jet Propulsion Laboratory 
4800 Oak Grove Dr. 
Pasadena, California 91109 
M.S.125-224 


Dr. Arnold J. Kelly 

Dept. of Mechanical & Aerospace Engineering 
Princeton University 

Princeton, New Jersey 08544 


Dr. Mitat A. Birkan 

AFOSR NA 

110 Duncan Avenue, Suite B115 
Bolling AFB, D.C. 20352-0001 


Dr. Ed Weiler 

Code SZB 

NASA Headquarters 
Washington, D.C. 20546-0001 


LCDR John Riggs 


10215 Centinella Dr. 
La Mesa, California 91941 


oh 








s¢ 





ee 


Se ee 


DUDLEY KNOX LIBRARY 
NAVA! POSTGRADUATE SCHOOL 
MONTEREY Un 93943-5101 


GAYLORD $ 


ALES Hird “ 
pd abgt ¢ 
te 







Supt Ey KNOX L = a 


AEN PMI 


ay a oe h D 
aa Roar Ret SS an ae q 
+ ie ake ‘ee 4 ohn Ar ay eat S 
ae ee aR eae eR a | 00307094 7 
e De $y A * q H te ba p b : 
é fue ee eee | i f c : 


atte res oe 


















© 
en 
Hf 


? 

ath eran 
i tad PEL cae ei 
re | ee tee eign gs = 
ce he * ' 
PANES et: 
' te Pre ees A 
ea Seed esa 


reg 















fd 
al hes ‘ot Kee 
Ai oe a 






oi 































7 re i OD ee ; 
. vs rid ew 
Pd Ma ete | nth 


‘ mens 
WO Cet 











Feed ey 
ia Me te 
seep ghee re 





Pa Pra’ 


dE ie P 







tho a | fei 
ar Laas) hi she 
Lyd PTUs 





you W 
. Re Hae Ae Res Hos ct) 
fete Te ee) Oe | 
Part it) Peat fee ie 
wee F +7 eee | oo Pa Ed ee 
Bios, Bhat iy a A tia Bhat U 
ane Gor ibe hid obs PLE eae W saahe a4 
? id Ae | fetes a Pa Naa 
rhs oe ret ie % eae ; 
5 


raat) . rie 
ust Soa yt E oir 4 <i om oe 

gaa Uehs Pata sf a be 
rreu Td Folin OF of fT 














qul %ea.ditoe rd a tear 
Nias pda id 
i Pe ery Ya 

a oe ahi 
oi sf te Stal 
rr 4 ris A i 
tet PH pe 
FN 4 

































Pret ae 
chal atid a 




















7 


? 








ore 






Pe 4c Ci i 
Vag ite tba he ; 
r ich eed ar 












AL STry st Ure ar 
aor ee ote te 













et er 


eer 









rs, Pogue et Pa ee 
iste Pa TACREM Ya nee aly 
7 Pag {1 




















Rh aie eee at) 
aaah 


pe wore) 
gi? Boe 






aH A ri 


Uri 
ee 
















Le 
aa rg Pe 
ahd P. Feheds 
Avetg * Me 
Py Ly atrsp° 
a Paar 
OE chet Late 


siete 
Ve ata 
Ary rol ba A ee 


oF pf CS 
reel EO AN Rae 
Le ae be 
“ee 































td 
ie " 














aes yee 
Pym 


rT ree 


J pees 

bat ry Ie IC Ur 
A aa) eed 
tat MEAT eer pa 5 




























ey a 









rae Cas ade | Pt ie a t « Et ite oer 
Fit gtels sapere’ ds ver Ba Aa ; Aree Apa ar cr Oren 
a ‘ Se » ee Lr 
Aer tet We 





i 
rd ai as oe 
fe ener e 












Bye ica 
he ¥ heir oe a 
ee ote 

re 


aot 
eae) 


Pa 
a a Pe a 2 










<r 5? ¥*_d te 
A gilinh ae 
rT 






















siete aS5 97 










. ee | 








Ca he aa ao Pr 4 a 
es {hf t ase) ‘a3, a ’ a aah 
> ty a PA ad Cr ae Sees ey i 






feb Pots 





Bi aes 


















































Pasay’ re ST es Pte 
PR ree S882 HE fhe a 
ry 
cB 
i 
ad oad ° 
paren pe 
be He ak Atl ] 
Pa Mt rv 
Ata Efe . or hea aa 
Brod Moa aM ele Bd AME NOR dey a 
5 ie iL See es or » A ro at 
rT a ed a Re ae Leal ie et i 
oer Ca eae dtd od By hat ea. i] Coheed - PT) *. 
ere a $2 A a 





2, SMF pres 6°"L9, Coe Rene” ¢. 8 
a bares tin web bgrzn Hh Mostys 
PUP TT bake Othe PY ert ee er ree ie 

U a ae ee aii re et eee ry se ogee 
eed a Lon) we" > ye 


iid. rk 
Pf iaege haa cl et 
fyew are aay 

CoE 5 eo aT 

























J “i o re ae .*% asisoverdy &° 
“Ar otek ae an tg ae i ee 
Ca i 











es | 

i hee £ he 1% Vim Bee ¢ 
oe 

ot Sd aun es be Lee) es | 

pee yt 2 ey 

ree a ere Mae a 








a 
dk ihe deed Hea 
Sal ded Ta td nad 
‘srobi dues! A Pa esos, 
RAM 4) deta ey Lat 

beat ee Ey 7 w 
Ce eet he Heresies 1) ade iW 
w betes Yeon) 

es rae ee 











































ed 7 
Sh es ff 
ort we Pie 











Pea nae my By a Pay a ce ee Sah er 
et a ye fo be al be et! a ee ca Ps o nt amare 
ey av uds tn “rs ts ee tery > eee oT iy 
Mieoers aA te typ es es Re ps Pei atc 
4 ‘ 






















Tire S 


ie a al PC ee er eee | 
Pays ue oot td 


veut iy toe aed a TE bapt-ocve FP 





Raa 
















Td ed Ae eat) ‘ ee ee ty het) Pe Py rt] te, & peas of Boar Cee } § 
9903" whe. Pe ee eee be eee ee ee ey fay aye 
TU ae eps rer % ae He . Jide WL nae whore, 





Pot ee ed * “o alfet 
eh a) adh 
























rr Ly © Bee 
Hee a : De ey 
~~ ze 
nd hee Led @ Cea 
a iors 












rid i yy 
Laie Ml Per or J + 
2 Pt omega om a * A tet eo 









iisere 
















it 
Te to at 


Pos ere 
weeSysta lp 274s 


tse ry Rrcrare ree) 
” Le ree erg 


i 
Ra blend 





Aaa nL 
pas® +403 440 


aa 
























¢ 
al ae ee er tit ry 
are ee) 




























hy tae rea a Th 
Bah Rants Ta 
fw ese Laren Ait rans 





“EF, Ss HRP vey oe ree ar are ar 

































































a ce Sae ye eat “waa * pie e's 
F ad - 
rar Au hae s oh Lia tA er hee * * Peete Tt | CL ° 
| a | bi hoe ices he aed Safar rrr aera tr) eer aria 
R k Panty i ror 7% ? . eee 
e ¢ >? Pope f) ry a 
are be tt og? errr 
oars feasts > 7 as i Aa Pitre sea 
ee eee ta Pe ae ~*~? sei ne * Prd} 
A 


F ves ni 

oS A aa a Pe 

aed el ee ho Soe 
7 am wsrgrg? ser it 

Sy) Radel SA betel ed i ad 

eu gory eee ats 

1 he bab Ag Tt ees tS td 

bie a4 & 
Wie i tes 2) PTE uh A ee ole 
a au ohh ald Le Ly 








od i 
meredae ye 













teil bd Be, iy 
a in pelrar 
opty cA 
bas 


















iad iy 
eee rh “Lr Fil 4 i rae i) pele ear 
Pred 




















































7 
a eae Heer re eat) ae LLY rer 
et Tt ho ee re A Be batt "3 ATR “ere des ipa t i : 
iss o bem oml ut hipts by ie oy CB. a ra F i Pee 
laa Te ee | a ry ee Pid ede S rt Be Be aerate 
haps Ped ree. ee ee eT) CP erie erat irs Pas BOERS atts Rt 
Td bel a be tk Se te wet kt te Saal cece ie He ere aa Pe Was fi ary id cs Pe ae 2 
ree ye de gees ea benatir Me eto y we eale gy a Pa ee ee eA pat sae Pia ma oy nt A 
Es iad Ca *e' zh “ wi PMS me | eae tT) Aa 7) ers a 
a ee | re “og tt: “apt B19 4 , 





ai 
Ly 





oD Cad ita 
ae 





“i Ae ha 












ry 
are IF doe Pots oy 
ee ed ‘to Le strait 
wt ' © eben ge Pot fhe 
“aR laa ey a, 
a rs we rae} Py 
fr ML to2 iy ai ee Ie fea 
id ead Pre 


7 






































" es ; 
Tad ee he he bs ats ri tomy ; jr 






























































pad fos Yes d ad M68 wae ti ry o> 
a St Bs Kote se t- a etre a i 
ih BR FLEES ee Bae er Eee erate Pies 
pA sO RYIM Tare sadgimdentar Opn e: S}eeHi eh yate-ve a wii le rveyed 1p reer TA 
o ak ve Ped oe re OT ieee Cae mad) mpeh lye (Co 3 4 Cae er be i 
LAAs he oes bb Bass bia Pane an ue ed Gan tak 2 Ae Ck We epreete of Cit vr: 
4 SSH itrD ere FL ee eyed Te Shae de sh 2 do o*test 8 ad Cree 
ao Penta hae Mem ue ie ute ub bh Lard Tn ae) 
"gigs Pr Pah ae a rhe ee) 
‘aes am Repay rl Pye eee 
peg LN a A 
id ue TR Br t Ue v 
x Sees hry . = A 
oT ht 7; * J 
*y Se) \ rhe to * hh 
Pe a | 


4 Mad Tedhale 
io A 4 ra Oe may 
Lee a er) i] 
a. o ley 
ae Ba . 


eh a 
Lf a 






¢ or reer 
oh, Por Ay i 

moh mE reese a! ie Perit me of) Porn td 

fas Pie Cobol. eo heer HLA ears ar oo 
Eas ry pera aoe? TT Pe erg ys ys ty trae 

Pele te nl ie PEt hota Oe id eed yea RH 





































, 
ea Wal . 





































al at yas r yererts Th J oi 
athe eine ea mR ree eae “aes 
; 5 Heres a iy A ar BU 
ho Tae amid, and x4 re yp RSet peat ays CES re ise Bie v cS = BY A ea Fame ce 
ar bate my us 7 ae 3 ere he om Nat } eee 5 Sand es 
>t Aes had rt ety ay to Lay " Ph eines ‘eee yi) i aniearerir 
Peal ae wer Eee U hited} a woes o Spa, OED a! hs ht e a aa reer ae we 
Cea es set ete ite Fy were RRP Mok COLE CC Ler mer oP bt. ee Sais Caley San ; ‘ae oa , 
eet eee See TERY ca es Pra He tale SO ey LaLa AA en Pees ey yar erik ae - ae 


‘ 
AUS th (ees a ter edgees ra Car i | 
aii mc a rh an dtey gs Ce at is PL a A 
7 Baa 


fry TY 
icoinua Ren 
WM i ag t 
he Aiea 












Lia 
Are Se 
he) ¥ gery eA 
eaNer ‘ 








= 














ar ta a phe 5 se . 7 P 
yee 3 meres an eT 

i> a one Pm oS 
Vel Yd pet eli} att} Bt La ro Bb AAS) rer apema hog 


